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CHAPTER I. 


INTRODUCTION 

Weather radars and imaging sensors on geostationary weather satellites are cur- 
rently the most widely used remote sensing tools for the short-term forecasting or now- 
casting of warm season convective storms and for warning of severe thunderstorm 
hazards. Zipser (1983) defines nowcasting as "the description of the state of the current 
weather and forecasts within the valid extrapolation range for each phenomenon which 
are based on intensive observations". The valid extrapolation range is further defined "as 
the period within which weather forecasts based upon observations and extrapolation are 
useful". The valid extrapolation period as well as the amount of lightning activity 
depend on the phenomena being described, geography, season of occurrence, instability 
of the atmosphere, and structure of the storm environment (Table 1). Dynamical 
forecast models with explicit physics are presently more applicable to greater time and 
space scales (Figure 1). 

In a discussion of the stages of nowcasting, Wilson and Carbone (1984) state "the 
first element of forecasting is simple extrapolation of event position and intensity. Pre- 
diction of completely new development or onset of dissipation of the existing event is a 
distinctly more ambitious nowcast objective". Forecasts of the future location and inten- 
sity of clouds, precipitation, lightning, or storm severity can be assessed by asking 
yes/no or how much. Did it rain at all ? Was there any lightning with that storm? Was 
there severe weather (flooding, hail, tornadoes, microbursts)? How much rain was 
forecast? How much lightning? Extrapolation forecasting is akin to conditional expecta- 
tion. What is the probability of rain at point P x in the next hour or two, given that it is 


Table 1. Typical Linear Extrapolation Time Scales for Various Weather Events 


Weather Event Time Scale for Linear Nonlinear Pre- Accompanying 

Extrapolation Validity dictive Capability Lightning 

(Nowcast) (Beyond Nowcast) Activity 


Downburst/Microburst 

-1 to a Few Minutes 

Tornado 

-1 to a Few Minutes 

Thunderstorm, Individual 

5-20 Minutes 

Severe Thunderstorm 

10 Minutes to 1 Hour 

Thunderstorm Organized 
on Mesoscale 

-1-2 Hours 

Flash-Flood Rainfall 

-1 to a Few Hours 

High Wind, Orographic 

-1 to a Few Hours 

Lake-Effect Snowstorms 

A Few Hours 

Heavy SnowAVinter 
Storm/Blizzard 

A Few Hours 

Frost/Freeze 

Hours 

Low Visibility 

A Few Hours 

Air-Pollution Episode 

Hours 

Wind 

Hours 

Precipitation 

Hours 

Hurricane 

Many Hours 

Frontal Passage 

Many Hours 


Very Limited 

Often Many Intracloud Flashes, 
Few Ground Flashes 

Very Limited 

Often Many Intracloud Flashes, 
Ground Flash Rates Variable 

Very Limited 

Variable Ratio of Intracloud to 
Ground Flashes 

Very Limited 

Typically Many Intracloud Flashes, 
Ground Flashes Variable 

Some 

Typically Many Intracloud and 
Ground Flashes 

Very Limited 

Varies From Many Ground Flashes 
to None 

Some 

- 

Very Limited 

Some 

Some 

Not Usually 

Some 

- 

Some 

- 

Some 

- 

Some 

- 

Some 

Variable 

Fair 

Variable 

Fair to Good 

Variable 


‘Adapted From Zipser (1983); Doswell (1986) 
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Time (h) 


Figure 1. Efficacy of different approaches to short-term forecasting over a range of 
time and space scales (after Dos well, 1986). 


raining at point P Q now? Any ability to determine the presence, increase, or decrease of 
lightning activity as a storm approaches a given location or facility will provide valuable 
information to the user of these data. 

With the recent advent of lightning detection and location technology, it is now 
possible to directly measure the existence and frequency of lightning activity in storms 
over large areas. The deployment of regional and national networks (Figure 2) using 
magnetic direction finding (Krider et al., 1976) and time-of-arrival (Bent and Lyons, 
1984) techniques to detect and locate cloud-to-ground lightning offers a new and 
decidedly different nowcasting data source for real-time multisensor data fusion. In the 
future total lightning rates (intracloud and cloud-to-ground) will also be observed in 
real-time by a lightning sensor in geostationary orbit (Christian et al., 1989). 

The observed lightning activity may be used in determining the existence, initia- 
tion, movement, dissipation, configuration, areal extent, intensity, and redevelopment of 
convective storms (Goodman et al., 1988a; Lewis, 1989). A recent evaluation of the 
operational use of lightning data by forecasters at the National Severe Storms Forecast 
Center (NSSFC) demonstrated great value in monitoring lightning activity for assessing 
the threat of existing storms and in issuing weather advisories; most frequently when 
storms were classified as strong (5-min update interval) and less frequently when storms 
were weak (1 5-min update interval). Furthermore, when forecasters were asked if 
lightning activity added knowledge about the general convective activity that could not 
be obtained from either satellite or radar data, a positive response was acknowledged for 
78% of 153 storm episodes considered strong and for 64% of 301 cases of storms con- 
sidered weak (as subjectively characterized by the forecasters). 

These preliminary results suggest that lightning activity and its association with a 
storm or complex of storms should be quantified and used as a source of nowcasting in- 
formation in knowledge-based and expert system/artificial intelligence/neural network 
algorithms being developed and tested (Browning and Collier, 1989; Roberts et al., 1990). 
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Figure 2. Demonstration national lightning network direction finder (DF) sites. 
Three sites are part of the MSFC network at Tullahoma, TN (DF 2 = Tl), 
Centerville, TN (DF 3 = Ce), and Barton, AL (DF 4 = MS). 

DF 1 (at MSFC) is not part of the national network. 
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The quantitative assessment of lightning activity offers a wide range of opportunities to 
develop algorithms to evaluate the synergism of these data as an adjunct to the weather 
radar, satellite, and conventional meteorological data (Watson et al., 1987; Goodman et 
al., 1988). The Advanced Weather Interactive Processing System (AWIPS) currently un- 
der development by the National Weather Service will be the first opportunity for all lo- 
cal forecast offices to integrate, process, and transmit high volume radar, satellite, upper 
air, and surface data. 

The Need for Real-Time Lightning Observations 

A recent survey of federal agency requirements provides a clear impetus for 
developing real-time techniques to monitor lightning hazards (MSI Services, 1986). Ex- 
amples include: 1) data requirements by the National Weather Service (NWS) during ac- 
tive thunderstorm periods to issue severe weather warnings; 2) timely information 
needed by the Federal Aviation Administration for flight safety, dispatch, and air traffic 
control/operations; 3) reliable tracking of lightning to support the safe and efficient 
operation of a variety of naval activities including operational and training flights, 
weapons and munitions handling, and aircraft and in-port ship refueling ; real-time dis- 
play of the location and direction of movement of cloud-to-ground lightning strikes 
within 300 km to support Air Force strategic and defensive activities including refueling 
operations, munitions handling, radar operations, computer operations, safety, and field 
exercises. 

Facility applications of lightning data include various test range activities such as 
those conducted by the Air Force and National Aeronautics and Space Administration 
(NASA) in support of the Space Shuttle and unmanned space vehicles, and by the 
Department of Energy in support of the Nevada test site where underground nuclear 
tests are conducted. Presently, NASA operates lightning detection systems at Marshall 


6 



Space Flight Center (MSFC), AL; Wallops Island, VA; and Kennedy Space Center, FL. 
Lightning data from a Navy operated network are also distributed to the Stennis Space 
Center, MS. The operational requirements for lightning data at test facilities are driven 
primarily by the concern over personal safety. However, improved lightning warnings at 
Kennedy Space Center also permit safe fueling operations during stormy weather periods 
at an estimated annual savings of one million dollars. Real-time users of the MSFC 
lightning network are shown in Table 2. 

Utilities presently use lightning data to design protection for power lines and dis- 
tribution systems, and to deploy repair crews during severe storms (Fischer and Krider, 
1982). In general, lightning warnings tend to be very conservative, with many opera- 
tional and training opportunities cancelled unnecessarily resulting in lost productivity. 
This brief summary of applications strongly suggests that any lightning sensitive tasks 
concerned with optimizing the safety and use of material and human resources can 
benefit from the currently available lightning detection and location technology. 

Review of Nowcasting and Extrapolation Forecasting Algorithms 

A nowcasting system consists of two main parts: 1) some type of characterization 
of the present weather situation and 2) a means (i.e., a model) to project the situation 
forward in time and space. Forecast methods using weather radar to extrapolate rainfall 
(an appropriate analog for lightning patterns), storm position, and intensity typically use 
some reflectivity (intensity) threshold to define the convective storms as either cells 
(clusters) or rain areas, correlate two or more successive observations to get a storm mo- 
tion vector, and extrapolate the intensity pattern some time into the future. 

The existing operational radar nowcasting systems extrapolate the characteristics 
and full intensity of the precipitation pattern without consideration for growth/decay of 
the rain intensity or its spatial distribution. The primary source of forecast error in a 
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Table 2. Real-Time Users of the MSFC Lightning Network 


U.S. Air Force, Arnold Engineering and Development Center, TN 
U.S. Army, Redstone Arsenal, AL Test and Engineering Directorate 
NASA, Marshall Space Flight Center, AL 

Earth Science and Applications Division 
Rocket Motor Test Areas 
Information Systems Office 

Safety, Reliability, Maintainability, and Quality Assurance Office 
Neutral Buoyancy Simulation Facility 
'Redstone Airfield Flight Operations 

WAFF 48, Huntsville, AL Television Station 

'National Weather Service Office, Huntsville, AL 

fState University of New York at Albany, National Lightning Network 


'Future Users 

tRaw Bearing Information Provided by Three MSFC Antennas 
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study of rain patterns associated with weather fronts in Great Britain was attributed to 
the development or decay of rain areas in 16 of 29 (55%) events (Browning et al., 1982). 
Storm growth/decay, mergers, splitting, and fragmentation also compromise the perfor- 
mance of peak reflectivity trackers (Crane, 1979; Rosenfeld, 1987), centroid trackers 
(Barclay and Wilk, 1970; Duda and Blackmer, 1972; Zittel, 1976; Bjerkaas and Forsyth, 
1980), and (pattern) correlation trackers (Austin and Bellon, 1974; Rinehart and Garvey, 
1978; Browning et al., 1982). 

Peak reflectivity trackers isolate and track local maxima in the reflectivity or 
precipitation field. Such techniques tend to overestimate the number of physically 
realistic storm cells, but do not miss the small, potentially severe storms that may not be 
identified by the other methods which depend on intensity thresholds to delineate 
storms. The centroid trackers use the 3-dimensional reflectivity weighted centroids 
(storm mass) to delineate , characterize, and follow storm movement. Correlation track- 
ers compare a field of reflectivity values at two successive times to get a motion vector 
for the entire system (more applicable to widespread light rain situations) or in localized 

sub-areas (applicable to individual thunderstorms). 

The Bjerkaas and Forsyth (1980) mass weighted centroid tracker has been imple- 
mented as a Next Generation Weather Radar (NEXRAD) system algorithm (NEXRAD, 
1985). NEXRAD is the Doppler radar system being jointly deployed across the United 
States and overseas by the National Weather Service, Federal Aviation Administration, 
and Air Weather Service to replace the aging WSR-57 and WSR-74 network radars now 
in service (Leone et al., 1989). These new radars offer numerous advantages over 
present weather radars in severe thunderstorm warning, rainfall estimation, and the 
detection of wind shears. 

In a study comparing the performance of different types of storm trackers, El- 
vander (1976) found the cross-correlation tracker performed best when presented only 
base-scan (lowest elevation level) reflectivity data and the centroid tracker superior for 
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volume scan data (the data acquisition mode for NEXRAD). A comparison between the 
Crane (1979) peak reflectivity tracker and centroid trackers showed general agreement of 
the storm motion vectors. Brasunas (1984) recommended using the correlation tracker on 
slowly moving widespread rain areas and the centroid tracker on the more convective 
storms. For large areas with multiple storm motion vectors. Browning and Collier (1989) 
suggested applying the correlation tracker to subareas within the confines of the larger 
system. 

Other sources of forecast error can be attributed to incorrect specification and 
delineation of the initial pattern (measurement errors), errors in estimating the initial 
pattern trajectory, and errors due to changes in the trajectory during the forecast period. 
Examples of this latter effect are storms that upon becoming severe tend to move to the 
right of the mean lower tropospheric wind (Newton and Fankhauser, 1964; 1975), and 
rainbands associated with tropical hurricanes and extra- tropical cyclones where the am- 
bient wind field imparts both a translational and rotational component to produce a cur- 
vilinear trajectory. 


Quantifying the Valid Extrapolation Ranee 

Linear extrapolation of the present trend is perhaps the easiest and most widely 
used method for nowcasting. However, the motion of the atmosphere and growth/decay 
of convective phenomena are examples of complex non-linear dynamical systems. The 
limitations of linear extrapolation with increasing time scale are readily apparent as the 
forecast error becomes large and then exceeds the threshold value for an acceptable or 
useful forecast (Figure 3). As noted earlier, the valid extrapolation period will depend 
on the process under study, which itself will have some mean lifetime. A non-linear 
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Figure 3. Possible extrapolation model fits to a non-linear process based on ^two . prior 
observations (OBI, OB2) and the current observation:a, linear extrapolation model, 
fS e^tetion range defined by accepiable error bounds; £, non-l.near 
model fit that is worse than linear fit; d, good non-linear model fit 


(after Doswell, 1986). 
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model may be inferior to a linear one if the nonlinear model is applied beyond the valid 
extrapolation range. In general, a non-linear model should provide the best forecast of a 
non-linear process such as the growth/decay of a thunderstorm. 

For large, approximately steady state weather systems containing widespread light 
to moderate rain showers, extrapolation forecasts can be very accurate for up to 6 h 
(Browning et al., 1982). Projecting the motion of intense thunderstorms that persist for 
several hours (supercells) is not too difficult, but predicting if, when, and where it 
might spawn a tornado is beyond present capabilities. The tornadic storm that struck 
Huntsville, AL on 15 November 1989 is one example of an isolated supercell storm that 
could be identified on radar in Mississippi 5 h before it reached the city and produced 
an F4 intensity tornado that killed 22 people. Yet, the tornado struck with almost no 
warning. Forecasts of convective weather phenomena at smaller space and time scales 
can be much more difficult since they change size, shape, and intensity more readily. 

Study Objectives 

The objectives of this study are twofold. First, the applicability and limitations 
of extrapolation techniques are explored relative to the problem of forecasting lightning 
(i.e., thunderstorm) activity at a future location and time using data acquired primarily 
from the NASA/MSFC ground strike lightning network installed, operated, and main- 
tained by the author since 1985. Radar echoes and the spatial distribution of the lightn- 
ing itself will be used to analyze, characterize, track, and examine extrapolative forecasts 
of storm position and the accompanying lightning activity. Second, the usefulness of 
physically-based non-linear models will also be examined for their applicability to the 
lightning forecast problem and to determine the valid extrapolation range for different 
weather scenarios. 
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In Chapter 2 the lightning and precipitation time series for storm systems over a 
wide range of space and time scales are examined. These results are used to develop a 
conceptual understanding of the thunderstorm life-cycle. The three parameter logistic 
model is offered as a candidate model of the thunderstorm life-cycle. Its properties and 
their relevance to the extrapolation nowcasting problem are examined. The success of 
an extrapolation forecast at time t+At is influenced by the quality of the data source 
(analyzed in Chapter 3), an accurate description of the present situation at time t 
(examined in Chapter 4), and correctly accounting for changes during the forecast period 
(investigated in Chapter 5). 

Chapter 3 discusses the process of finding the most accurate locations of the 
cloud-to-ground lightning discharges. This process involves removing systematic errors 
from the data and the application of an optimization technique to locate the most prob- 
able position of each lightning discharge. A novel approach using isolated radar echoes 
to constrain the error correction and optimization procedures to remove systematic errors 

is described. 

Chapter 4 presents a pattern recognition scheme that is used to generate initial 
seeds or "first guess" fields for clustering the discrete lightning discharges into storm 
cells. The clustering process is critically dependent on the prior accuracy of the lightn- 
ing location estimates (Chapter 3) and the generation of subsequent storm life-cycle time 
series (Chapter 5) relies on the cluster analysis procedure assigning the correct number 
of lightning discharges (objects) to the proper storms (groups). The advantages and 
limitations of different clustering strategies for storm identification and tracking are ex- 
amined. Storm identification with lightning data alone is compared to storm identifica- 
tion with radar alone, and some synergies for sensor fusion are explored. 
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In Chapter 5 the logistic growth model is utilized to examine the storm life-cycle 
over a wide range of space and time scales and address the potential of non-linear 
regression models to improve upon short-term extrapolation forecasts. A physical inter- 
pretation of the logistic model parameters and the resulting implications for determining 
the valid extrapolation range is considered. 

Chapter 6 summarizes the chief results of this research and discusses how these 
results move the state of knowledge forward. Chapter 7 concludes the discussion and 
offers future areas for additional research. 

Appendix A provides the mathematical formulation for the optimization algo- 
rithm (not currently available in the open literature) employed in Chapter 3 to find the 
most probable flash locations. Appendix B contains the FORTRAN-77 source code used 
to convert the raw data into geophysical quantities, correct the systematic errors in the 
data, and compute the optimal flash location. Lastly, Appendix C contains the 

FORTRAN-77 source code for the cluster analysis. 
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CHAPTER II. 


A CONCEPTUAL MODEL OF THE THUNDERSTORM LIFE-CYCLE 

A conceptual understanding or model of the thunderstorm life-cycle is needed 
for addressing the utility and limitations of extrapolation forecasts of non-steady-state 
weather phenomena such as thunderstorms and their accompanying lightning activity. 
The following discussion examines the lightning and precipitation time series for storm 
systems ranging in size from an isolated airmass storm to a large mesoscale storm com- 
plex encompassing an area of nearly 12,000 km 2 . 


Air Mass Thunderstorms 


The relationships between the early electrical development of a thunderstorm cell 
and the vertical development of the radar echo, precipitation, and cloud top are depicted 
in Figure 4. A significant fraction of the total lightning (50-95%) occurs between 
regions of opposite charge within the cloud (intracloud lightning) without ever reaching 
the earth (cloud-to-ground lightning). The initial discharge will almost always be 
intracloud. Typically this discharge occurs 5-10 min after initial electrification, which 
itself begins 5-10 min after the detection of a 35-40 dBZ radar echo aloft. Based on 
lightning and radar data collected in three different climatic environments (New Mexico, 
Alabama, and Florida), Buechler and Goodman (1990) find that the time lag from the 
reflectivity exceeding 40 dBZ at the -10°C level (the height of the main negative charge 
region in the thunderstorm central dipole) to the first intracloud discharge ranges from 
4-33 min. The time lag is related to the rate of vertical development of the cloud. On 
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Figure 4. Composite chart of thunderstorm development observed by radar, radiosonde 
and electrical sensors (after Workman and Reynolds 1949). 
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average, the first cloud-to-ground discharge will be of negative polarity (it lowers nega- 
tive charge to ground) and will occur 15-20 min after the initial radar echo is observed 
in a vertically growing cloud. The stage is now set for the active lightning phase of the 
thunderstorm life-cycle. 

As the cell continues to develop, the active electrical phase may have a duration 
lasting from a few minutes to many hours. The total amount of lightning and peak flash 
rates of a storm are a non-linear function of its height, size, mass, duration, and en- 
vironment (Shackford, 1960; Livingston and Krider, 1978; Williams, 1985; Cherna and 
Stansbury, 1986; Goodman and MacGorman, 1986; Goodman et al., 1988b). 

Figure 5 shows the relationship between lightning occurrence and precipitation in 
a small airmass thunderstorm 26 km 2 in area (>18 dBZ) observed near Huntsville, AL 
by the NCAR CP2 radar on 20 July 1986. The storm produced a strong microburst with 
a velocity differential of 30 m s' 1 . Total lightning (intracloud and cloud-to-ground) ac- 
tivity was measured by an instrumented mobile laboratory operated by the National 
Severe Storms Laboratory (Rust, 1989). The mobile laboratory was also used for ground 
truthing the lightning strike network, discussed in greater detail in Chapter 3. 

The 20 July case represents a Byers and Braham (1949) type airmass thun- 
derstorm. The mobile laboratory was situated under the storm throughout its 45 min 
life-cycle and recorded 110 intracloud Hashes and 6 cloud- to-ground flashes, all 6 of 
which were detected by the ground strike network. The first intracloud discharge was 
observed about 4 min after hail was initially indicated by radar, during a period of rapid 
vertical development as the cloud top neared its maximum height of 14 km. The first 
ground discharge occurred 5 min later when the maximum reflectivity core descended to 
5.5 km and a weak outflow was detected by the radar. 

Storm rain rates are computed from empirical Z-R relations developed by Mar- 
shall and Palmer (1948), Jones (1956), and Seliga et al. (1986). The storm ramflux (kg 
s' 1 ), mass (kg), and vertically integrated liquid water content or VIL (kg m 2 ) (Greene 
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Figure 5. Lightning and precipitation history of an airmass thunderstorm that produced 
a microburst on 20 July l^o6 at Huntsville, AL: above left , total flash rate time series 
where N. c and N CQ represent the number of intracloud and cloud-to-ground flashes 
produced by the parent storm; right , time- height cross-section of cloud top infrared 
temperature (°C), 0 dBZ and 30 dBZ echo contours and maximum reflectivity; lower 
ls£t, vertically integrated liquid water content (VIL), storm mass, rain flux and echo 
volume; right, peak reflectivity (Z) and storm average rain rates (R) as indicated 
(after Goodman et al., 1988). 
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and Clark, 1972) are calculated using a 30 dBZ threshold and the Jones (1956) relations 
Z=486R X 37 and M=0.052R° 97 , where Z (mm 6 m' 3 ), or Z H , is the CP2 reflectivity at 
horizontal polarization, R (mm h' 1 ) is the rain rate and M (g m 3 ) is the liquid water 
content. 

The peak total flash rate of 23 min' 1 was reached another 3-4 min after the ini- 
tial lightning, 6 min prior to the maximum microburst outflows, and in conjunction with 
the peak in vertically integrated liquid water content (5.3 x 10 s kg m 2 ), echo volume 
(1.9 x 10 11 m 3 ), and storm mass (3.3 x 10 8 km). The rainflux (1.5 x 10 5 kg s x ) and 
storm averaged rain rates (18.2-28.2 mm h' 1 ) reached their maximum values in associa- 
tion with a visual confirmation of pea-sized hail mixed with heavy rain about 2 min 
after the peak flash rate. 

Rapid-scan (5-min interval) satellite imagery from the GOES-E geostationary 
weather satellite were collected during the storm life-cycle. The infrared temperature of 
the cloud continued to show cooling (which could be misinterpreted as continued vertical 
development) even as the radar echo top and lightning rates decrease (indicating storm 
collapse). The misleading satellite signature is due to the small size of the storm and the 
(effectively 4 km x 8 km) field of view of the infrared radiometer (which also "sees" the 
earth’s surface). The radiometer field of view is underfilled at 1900 UTC. The field of 
view becomes more fully filled as the anvil expands, thereby sensing a decreasing cloud 
top temperature. However, the abrupt decrease in the total flash rate indicates storm 
collapse and thus serves as a microburst precursor. Yet, no such signature exists in the 
cloud-to-ground lightning evolution due to the small number of events (samples). 
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A Multi-cellular Thunderstorm 


Figure 6 shows the cloud-to-ground lightning and convective rainflux calculated 
every 10 min from the WSR-57 radar at Nashville, TN for a multi-cellular storm ob- 
served on 25 July 1986 over a period of 90 min. The convective rain area is simply 
defined here as the precipitating area within the 30 dBZ reflectivity contour. The 
lightning and rainflux are in-phase and are fairly well correlated (r=0.77). However, as 
the storm decays the lighter rainfall area contributes more to the total precipitation such 
that there is more rainfall per flash during storm decay than during storm growth. 

Convective Storm Complexes 

Figure 7 shows cloud-to-ground lightning and rainflux during a 7 h period of 
observation on 13 July 1986 of a mesoscale convective system in the Tennessee Valley 
that develops an extensive trailing stratiform rain region in the latter part of its life- 
cycle. The lightning data recording was briefly interrupted for a tape change at 2240 
UTC and continued at 2248 UTC, but the latter period is not shown here. Not long 
after 2300 UTC the entire storm system could not be sampled adequately from the 
Nashville radar as the storm moved out of range to the south. This case again shows ex- 
cellent agreement (r=0.96) between the lightning and convective rainflux time histories. 

Figure 8 presents a summary of the cloud -to-ground lightning and rainfall time 
histories of mesoscale convective complexes (MCCs) in the Central United States studied 
by Goodman and MacGorman (1986) and McAnelly and Cotton (1986). Such storm sys- 
tems are readily identified by their persistence and extensive cold cloud shields in in- 
frared satellite imagery (Figure 9). The typical precipitating lifetimes of MCCs are on 
the order of 12 h with spatial extents of a few hundred kilometers. Much of the warm 
season rainfall in the Northern Hemisphere (up to 70% in the major crop growing 
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Figure 6. Lightning-rainfall relationships for a small multi-cellular 
storm on 25 July 1986. 


21 




Flashes (10 min)' 1 



7kM(UTq 


(10 min) H 


^UJuIy 1 ^ 1 9 8 8 h 6 tnin8 ' rainfaI1 relationships for a mesoscale convective system 


22 





Figure 8. Hourly cloud-to-ground lightning rates and rainfall as a function of 
Mesoscale Convective Complex (MCC) life-cycle (adapted from Goodman and 
MacGorman, 1986; McAnelly and Cotton, 1986). 
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Figure 9. NSSL lightning network shown in satellite projection with the infrared 
cloud top image of a MCC during its mature phase (maximum cloud shield extent). 
The 4-DF deployment during 1983 and the 350 km range ring are denoted by crosses 
(after Goodman and MacGorman, 1986). 
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regions of the U.S.), extensive flooding, and severe weather is a result of these organized 
mesoscale circulations (Fritsch et al„ 1986). Up to 25% of the entire annual lightning 
strikes at a given site can be accounted for by the passage of just one MCC (Goodman 
and MacGorman, 1986). Due to their meteorological and economic significance, any 
nowcast/forecast skill that can be demonstrated for mesoscale storm systems is a 
worthwhile endeavor. 

The cloud-to-ground lightning activity increases and decreases exponentially over 
the life-cycle of MCCs. Based on earlier studies and the preceding analysis, this 
relationship appears to be generally valid for isolated storms, multi-cellular storms, and 
the ensemble convection embedded in organized mesoscale convective weather systems. 
The high correlations (>0.9) between lightning and rainflux extends over three orders of 
magnitude from 10 1 - 10 4 km 2 . Earlier scaling studies indicate that precipitating cloud 
dimensions are self-similar over 5 orders of magnitude (Lovejoy, 1982). Clearly, the 
non-linear physical interactions that produce the microphysical and dynamical properties 
of clouds are also relevant to their electrification. 

Positive Polarity Cloud-to-Ground Discharges 

Positive polarity cloud-to-ground discharges are often observed during the dis- 
sipation phase of the storm (Krehbiel, 1986). In addition, positive polarity flashes fre- 
quently occur 1) from thunderstorm anvils and storms which become severe and produce 
mesocyclones, large hail, or tornadoes (Rust, 1986); 2) in association with long-lived wet 
microbursts in low shear environments (Buechler et al., 1988); 3) in the trailing 

stratiform rain region of mesoscale weather systems (Rutledge and MacGorman, 1988); 
and 4) in the northern section of mesoscale systems, aligned with the geostrophic wind 
and downwind from the most vigorous convection, which is dominated by negative 
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polarity ground discharges (Orville et al., 1988; Engholm et al., 1990). The mesoscale 
system that led to the Huntsville tornado produced positive polarity discharges in each 
category listed above during some portion of its life-cycle. 


Figure 10 summarizes these lightning observations into a conceptual model of the 
growth and decay of a typical thunderstorm and its associated total lightning activity. 
The temporal evolution of the lightning activity is in-phase with the development of the 
storm updraft and is strongly coupled to the life-cycle of the thunderstorm described 
above and in earlier studies by Byers and Braham (1949), Workman and Reynolds (1949), 
and others. These results show the electrical development of the cloud is intimately con- 
nected to its dynamical and microphysical development. Laboratory measurements by 
Jayaratne et al. (1983) suggest that the charge transferred per collision is a complicated 
function of particle size and type, cloud liquid water content, temperature, and even 
chemical composition. A possible inference from these observations is that the greater 
the production rate of precipitation and ice particles in a cloud, the greater the charging 
rate of the storm. This is partially supported by the growing success of numerical 
models in simulating the initial electrification of small thunderstorms (e.g., Ziegler et al., 
1986; Helsdon and Farley, 1987). Multi-cellular storms will exhibit impulsive updraft, 
downdraft, precipitation, and flash rate growth and decay. Thus, the time rate-of- 
change of flash rates also provides a signature of the growth and decay of the thun- 
derstorm. 
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RELATIVE VOLUME FLASH RATE (MIN' 1 ) 

IN CLOUD VERTICAL VELOCITY 


CONCEPTUAL EVOLUTION OF THE 

ELECTRICAL DYNAMICAL AND MICROPHYSICAL PROCESSES INSIDE THUNDERSTORMS 



STORM DURATION (MIN) 

Figure 10. Conceptual model of the temporal evolution of cloud electrical, 
kinematic, and microphysical development. 
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The Logistic Growth Model 


Hald (1952) states that "the function used to represent the relationship between 
variables should as far as possible be chosen on the basis of professional knowledge 
about the problem under discussion and the reasons advanced for this choice are of fun- 
damental importance as regards confidence in extrapolations". The storm system size and 
precipitation particle population have been shown to be important factors in maintaining 
the charging/discharging process. One can attempt to characterize this process by simple 
first order differential equations which have been applied to population dynamics to 
describe the phenomena of growth and decay (Hald, 1952; Bard, 1974; Boyce and 
DiPrima, 1977; Haberman, 1977). 

The type of model needed depends on the the type of growth that occurs. These 
types of models are mechanistic in nature, rather than empirical. Mechanistic models 
are derived from assumptions on the type of growth, and these assumptions can be rep- 
resented by differential or difference equations (Draper and Smith, 1981). Empirical 
models are chosen to approximate the unknown mechanistic models. One likely can- 
didate mechanistic growth equation is the "logistic" or sigmoid curve (Figure 11). The 
logistic curve has frequently been used to describe the growth rates of populations (cells, 
human and animal populations, chemical kinetics, telephone subscribers, and business 
transactions). 

Let t denote time or the magnitude of a growth factor which influences the size 
y of the phenomenon observed, then dy/dt denotes the rate of growth per unit time. 
Let the process be characterized by the general equation 

dy 

— = f(t,y) (2.1) 

dt 

where the growth rate depends on both time and the size of the population. 
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Consider the following three special cases where 


— = f(y)g(t). (2.2) 

dt 

Letting f(y)=l, y, and y(a-y) gives 

dy 

-7- - 8(0 (2.3) 

dt 


dy 

— = yg(t) (2.4) 

dt 


dy 

— = y(«-y)g(t), (0<y<a). (2.5) 

at 


In (2.3) the growth rate y depends on time, but not on the size reached. In (2.4) the 

growth rate is proportional to the size reached and to a function of time. In (2.5) the 

growth rate is proportional to both the size reached and the remaining size, as well as a 

function of time. The latter case is the one of most interest. Now, write (2.2) as 

dy 

— = g(t)dt. (2.6) 

f(y) 

Introducing the "logarithmic differential coefficient" (Hald, 1952, p.659) 


dln(y) 1 dy 

dt y dt 


(2.7) 


and letting f(y) = y(a-y) and g(t)=/9 gives a relation where the growth rate is propor- 
tional to the size of the population and remaining size (where a denotes the growth is 
limited to some maximum amount, i.e., the value that y approaches as t increases) as 
well as a function of time. The growth rate relative to its present size, 1/y (dy/dt), 
decreases linearly as y increases. The resulting solution of the differential equation is 
the logistic function 


a 

(I+/?e kt ) ‘ 


( 2 . 8 ) 
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This model is but one of many possible exponential growth models having similar forms 
(Hald, 1952; Williams, 1959; Richards, 1959; Draper and Smith, 1981; Bard, 1974; 
Haberman, 1977; and Myers, 1986). At t=0, the starting growth rate is a/(l+/9). Also, as 
t — » oo, y -» a. 

The slope of the logistic curve is positive and the second derivative (Draper and 
Smith, 1981; Prof. Don Ryan, personal communication) 

d 2 y k 2 

- 5 - = (— )y(<*-yX<*-2y) (2.9) 

dt 2 a 2 

has inflection points at y=0, a, and a/2. At the point of inflection y = a/2, substitution 
in Eq. (2.8) gives the time of inflection as tj=(ln/9)/k. We note that the curve is sym- 
metric about this point, i.e., the system decays or diminishes at the same rate at which it 
grows. Thus, for nowcasting purposes one might first compute the rate at which the 
lightning activity increases (i.e., a growth rate) and the time required for a storm to 
reach its peak discharge rate (the point of inflection). Based on symmetry, one would 
then predict the storm to decay at this same rate and reach the end of its life-cycle in 
the same number of time steps needed to produce the initial 50% of the total lightning. 

In order to test and evaluate this model, one must develop a methodology for as- 
sociating the discrete lightning events with their parent thunderstorms. This process is 
addressed in Chapter 4. Next, generate a time series at uniform sampling intervals. 
During each successive sampling period the parent storm must be tracked with time and 
correlated with its past position. The extrapolation forecast results are presented in 
Chapter 5. However, the pattern recognition process is strongly dependent on the 
quality and limitations of the lightning strike data which is described next in Chapter 3. 
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CHAPTER III. 


OPTIMIZATION METHODS FOR LOCATING LIGHTNING FLASHES USING 
MAGNETIC DIRECTION FINDING NETWORKS 

Introduction 

Magnetic direction finding (DF) networks for locating lightning strikes to ground 
require that two or more receivers detect the characteristic radio signal produced by 
return strokes (Krider et al., 1976). Once the signal is detected, an estimate of the most 
probable flash position, sometimes referred to as the best point estimate (BPE), and a 
confidence region can be constructed. The spatial distribution (or clustering) of the 
lightning flashes, however, is a function of the dimension and vigor of the storm, the 
orientation of the lightning channel and hence its radiation field, and the errors (both 
random and systematic) associated with the technique. 

The systematic errors due to DF site effects are a major source of network 
degradation (Ross and Horner, 1952; Horner, 1954; Gething, 1978). When one or more 
of the network DFs do not detect the flash, the location estimate must be determined 
from a less favorable geometry (e.g., a flash along the baseline of two DFs, or a flash 
more distant from one site than another site in a more optimal geometry). The 
reliability of a fix (i.e., position estimate) can generally be maximized by using only the 
two closest stations to the target (Stansfield, 1947). However, near the baseline of the 
two DFs it is better to use a more distant receiver in a more favorable geometry (which 
will usually produce a smaller confidence ellipse). This is the basis for the real-time al- 
gorithm implemented in the earlier versions of the lightning DF networks manufactured 
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by Lightning, Location, and Protection (LLP), Inc. (Krider et al., 1980). This algorithm 
chooses the DF pair having the greatest signal strengths (presumably the two closest DFs 
to the flash). If the flash is near the baseline of the DF pair, a solution can be com- 
puted with the DF having the next strongest signal strength. An algorithm called 
"multiple correlation optimization" now replaces the simple 2 DF technique described 
above when three or more DFs detect a flash (LLP, Inc., 1988). This algorithm basically 
performs a least squares minimization between the most probable flash position and the 
sum of the bearing errors. 

This latter algorithm, first introduced by Hiscox et al. (1984) and a more recent 
eigen-vector algorithm introduced by Orville (1987) are very similar to a technique first 
proposed more than ten years earlier by Wangsness (1973). All three algorithms attempt 
to minimize the same objective function although different methods are used to reject or 
flag "wild" bearings. Hiscox et al. (1984) also proposed the use of properly normalized 
signal amplitudes as an additional weighting factor (or constraint) to determine the op- 
timal location. The improvement in solution accuracy by this latter method is a function 
of network geometry and the lightning location relative to that geometry. More 
recently, stochastic optimization techniques known as simulated annealing (Kirkpatrick 
et al., 1983; Szu and Harley, 1987a; 1987b) have been successfully applied to the general 
multiple DF/multiple bearing problem. 

This chapter describes the application of an eigen- vector algorithm (called FFIX) 
which is based on the technique described by Wangsness (1973). FFIX was developed on 
or before January, 1973 (but apparently never published) for finding radio transmitter 
locations anywhere on earth from multiple DF bearings. The technique characterizes the 
measured bearings by bearing planes that pass through the center of the earth and by 
their unit normal vectors. The BPE is determined from the vector from the center of 
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the earth that minimizes the weighted sum of squares of its inner products with the nor- 
mal vectors. A BPE and confidence ellipse are computed from the bearing data in terms 
of the eigenvalues and eigenvectors of a 3 x 3 matrix. 

This paper offers the first ever adaptation of the FFTX algorithm to the lightning 
location problem. The author is indebted to Dr. R. Johnson of the Southwest Research 
Institute and his sponsors at the Department of Defense for providing some documenta- 
tion and the source code in 1985. The mathematical formulation for FFIX is provided 
in Appendix A. Appendix B contains the FORTRAN-77 source code for the FFIX sub- 
routine. 


Overview of the FFIX Algorithm 

Individual DF bearings are corrected for systematic errors and correlated in time 
for each lightning discharge before being submitted to FFIX to determine the most 
probable lightning ground strike point (Figure 12). The input data consist of n station 
bearings and their respective standard deviations (random plus any remaining systematic 
errors). The root mean square (RMS) bearing error for each LLP DF is about 1°. Pre- 
vious attempts to iteratively remove the systematic error have not wholly eliminated 
them (Mach et al., 1986; Schutte et al., 1987). These techniques reduce the total bearing 
error (i.e, the standard deviation) to a value approaching 2° at best. It will be shown 
below that the standard deviation impacts both the BPE calculation and the confidence 
ellipse. 

Sines of the bearing errors are assumed to be independent normally distributed 
random variables with zero means, but some bearings may be "wild". "Wild" bearings are 
rejected by a process where all combinations of bearings are exhaustively evaluated until 
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Figure 12. FFIX algorithm flowchart. 
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a consistent subset of bearings is obtained. A set of n bearings is "consistent" if the as- 
sociated sum of squared bearing errors from the best point estimate is less than the 80% 
value of X 2 n _ 2 - 

The BPE and confidence ellipse are computed from the largest "consistent" subset 
of bearings. The algorithm employs two approaches for rejecting the "wild" bearings. 
The first approach, called the "exhaustive method", is invoked when ten or fewer bear- 
ings are submitted for a BPE. If the submitted bearings lack sufficient consistency to 
form a BPE, then all subsets of n bearings are taken (n-j) at a time, where j= 1,2,3,..., 
until an acceptable solution is identified or a lower limit on the number of bearings is 
reached, in which case there is no solution. The lower limit is the greater of n/2 or 3. 
Thus, the largest subset of consistent bearings forms the BPE. A "sequential method" is 
employed when more than ten bearings are submitted for a BPE. If all n bearings fail 
the consistency test, then the bearing that is more "inconsistent" with the bearing set is 
rejected (i.e., the wild bearing). The remaining set of (n-j) bearings is examined as 
before. In this way the most "wild" bearings are rejected sequentially. 

In the event that the chi-square consistency test fails, the "best solution" is 
chosen as the intersection of the bearing pair having the minimum semimajor axis in its 
confidence ellipse. This argument assumes the best network geometry also gives the best 
solution, all other factors being equal. This approach is supported by the earlier work of 
Stansfield (1947). This final iteration is necessary in a small direction finder network 
such as that run by MSFC (4 DFs) because nearly 50% of all flashes are seen by only 2 
DFs. The probability of detection for the network as a whole would be significantly 
poorer without this last iteration process. 
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Proof of Concept 


BPE Calculation for a 4-DF Network 

The Marshall Space Flight Center (MSFC) 4-DF network has been in operation in 
the Tennessee Valley (southern TN and northern AL) since 1984. Figure 13 shows the 
deployment of the DF stations and other ground-based remote sensing systems such as 
radars and rawinsondes. The additional systems were operational during June and July 
1986 in support of the Cooperative Huntsville Meteorological Experiment (COHMEX) 
multi-agency field program (Dodge, et al., 1986). The radar and rawinsonde data are 
used in this study to determine the location and characteristics of thunderstorms as well 
as the structure (vertical profiles of temperature, humidity, and winds) of the storm en- 
vironment. 

At the present time FFIX is used only in post analysis and follows the exhaustive 
rejection path shown in Figure 12. An example of the information computed for each 
flash is given in Table 3. The flashes are stored in rows and columns which indicate the 
hour of occurrence and flash sequence number. The file also contains the Julian day 
(DAY); hour (UTC) of occurrence (TIME); number of flashes in the selected interval 
(CMAX); flash time in hours, minutes, and seconds (HMS); latitude (LAT) and longitude 
(LON) of the flash; an estimate of the first stroke peak current in kAmps (KA); the 
maximum normalized signal strength (NSTR); the number of return strokes (RS); the 
semimajor axis (SMA) in km; semiminor axis (SMI) in km; orientation (ORI) in degrees; 
and area (AREA) of the error ellipse in km 2 ; flash polarity (FLAG= V for negative and 
V for positive); and time of the flash in hundredths after the last second (MS). The 
flash polarity is also indicated by the sign of the KA and NSTR fields. 
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LONGITUDE (DEGREES) 


Figure 13. Deployment of the lightning, radar and rawinsonde network for 
the Cooperative Huntsville Meteorological Experiment (COHMEX) conducted 
near Huntsville, AL in 1986. 
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Table 3. Lightning Discharge Data Base 
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Figures 14 and 15 show an expanded view of a radar echo (as seen from the CP-2 
radar) and the accompanying lightning strikes (shown as dots enclosed by the 50% error 
ellipse) for a small storm in Tennessee on 13 July 1986. The ellipse is generated from 
the type of information contained in Table 3. The lightning centroid is less than 5 km 
from the main echo. Figure 15 is produced using a bearing standard deviation of 2°. 
Figure 16 shows the corrected locations, but with a 1.5° bearing standard deviation, su- 
perimposed on the radar reflectivity. The spatial dispersion of strikes is on the order of 
10 km and is well correlated with the radar echo. By reducing the bearing standard 
deviation from 2° to 1.5°, the semimajor axes of the 50% ellipse for flashes 6 and 7 
(Table 3) are reduced (i.e., improved by) 40% and 24%, respectively. 

Earlier studies of the natural distribution of lightning strikes produced by iso- 
lated storms show that the lightning clusters are approximately 10 km in diameter 
(Feteris, 1952; Hatakeyama, 1958; Krider, 1988; Goodman et al., 1988). In contrast, 
however, the spatial distribution of the lightning in the trailing stratiform rain region 
behind summertime squall lines and in wintertime regimes can be very sparse and 
widespread (Rutledge and MacGorman, 1988; Engholm et al., 1990). 

BPE Calculation for a 6-DF Network 

Figure 17 shows an example of a flash detected by the 6 (now 7) DF lightning 
network operated by the National Severe Storms Laboratory (NSSL) in Norman, OK 
(Rutledge and MacGorman, 1988). This example shows how the consistency criterion is 
used to reject "wild" bearings. The BPE labeled FIX #1 (SMA=9.2 km) uses all six DFs 
as input, but the BPE uses only five bearings (DF 6 is rejected). The BPE labeled FIX 
#2 (SMA=32.1 km) uses only the four DFs in Oklahoma as input, using all four. FIX 
#2 using only the Oklahoma DFs is located 28 km northwest of FIX # 1 . The other 
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Figure 14. Lightning during the period 2315-2325 UTC on 13 July 1986 
superimposed onto the CP-2 radar echo prior to ground truth corrections: 
= negative polarity discharges; + = positive polarity discharges. Radar 
reflectivity contoured at 18 dBZ and 40 dBZ. 
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Figure 15. The location and 50% error ellipses (2° bearing standard deviation) 
for each discharge after radar ground truth corrections. 
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Figure 16. Overlay of corrected lightning location estimates with the radar echo. 
Same as in Figure 15 but with 1.5° bearing standard deviation. 


43 




Figure 17. FFIX solutions for cloud-to-ground flash detected in Kansas by 
the NSSL DF network. 
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intersection points where bearing pairs intersect represent false or ghost targets. Figure 
18 also shows the corresponding radar echo distribution and the location of FIX #1 
(indicated by the large cross in southwestern Kansas). 

Table 4 gives the error ellipse characteristics for FIXs #1 and #2, and for each 
bearing pair. When a BPE cannot be computed for a. = 1.5°, one can increase cr by 0.5 
increments until a solution is acquired. In this way, the standard deviation is treated 
more as a tolerance factor, rather than as an absolute. 

The results of this test show that FIX #1 has the smallest error ellipse of all sub- 
mitted DF combinations. For all paired-bearing combinations the DF (1,7) BPE has the 
smallest semi-major axis, but the DF (4,7) BPE has the smallest ellipse area. Although a 
solution is possible near the baseline of DFs (2,3), the error ellipse is quite substantial 
because of the poor geometry, reflected in the large value of tr.. The three greatest sig- 
nal strengths are reported at DFs 4, 3 and 6, in that order. An algorithm that uses the 
bearing pair having the greatest signal strengths would form a solution using the DF 
(3,4) combination, but the resulting semi-major axis for the ellipse is 201.9 km and the 
solution is displaced 14 km from the FIX #1 BPE. A solution could not be obtained 
until a. was increased to 5°, again indicating a poor geometry. However, a solution 
formed with the DF (4,6) combination reduces the semi-major axis to 23.3 km. Yet, 
any solution with DF 6 must be treated with caution, since it is inconsistent with the 
other bearings. One should ask whether the DF 6 bearing deviation is due to site errors 
or if it is just a "wild" bearing. If this particular DF 6 bearing angle is inconsistent with 
other flash BPEs in the same direction, then the deviation is due to systematic or site er- 
rors. However, if it is just a chance occurrence, say 9 of 10 bearings at this angle are 
consistent with the other DFs, then it was most likely a "wild" bearing. Following the 
completion of the aforementioned analysis a large systematic error at DF 6 was inde- 
pendently confirmed and the station moved to a new location in northern Kansas. 
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Figure 18. Distribution of radar echoes and lightning activity corresponding 
to the time of the FFIX solution in Kansas. Large cross indicates solution #1. 
Radar contours at 18, 30, 40, 45, 50 and 55 dBZ. 
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Table 4. Error Ellipse Characteristics for 6-DF Network 


DFs 
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a 
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LONG 

(deg) 

SMA 

(km) 

SMI 

(km) 
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37.557 
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1904.0 
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System a ti c E rror Correcti ons 


The FFIX algorithm has also been used in this study to correct for the systematic 
errors arising from site effects (Horner, 1954; Gething, 1978). Site errors cause bearing 
errors that are themselves a function of direction. The 12° bearing deviation from the 
BPE at DF 6 in the previous example was due to unresolved site errors. Site errors are 
one of the chief limitations to achieving the optimal location accuracy, and perhaps the 
most difficult problem degrading network performance. 

Previous attempts to correct for the systematic errors associated with lightning 
direction finding systems employed some type of optimization procedure that minimized 
the difference between the observed bearings and the "true" target location. The "true" 
target location can be determined by visual ground-truth (Mach et al., 1986) or by as- 
suming that one or more of the direction finder bearings is correct (Hiscox, 1984; Or- 
ville, 1987). This latter method will give self-consistent solutions (as applied to the DF 
6 bearing deviation described above), but spatial bias effects may still be present (e.g., 
lightning clusters offset from radar echoes). Rocket triggered lightning strikes at Ken- 
nedy Space Center, FL have also been used to provide ground-truth, but this is only ap- 
plicable to a single bearing line from any direction finder station. In practice, it is very 
difficult to obtain visual ground-truth at a large number of discrete bearings. Schutte et 
al. (1987) reversed the role of transmitter and receiver by radiating a 1 MHz signal suc- 
cessively through each loop of the direction finder antenna whose amplitude could then 
be measured by a radio receiver at a number of bearings from the site. 

The ultimate litmus test of any of these methods (and the practical usefulness of 
these data over a large domain) should be determined by the degree of spatial correlation 
between clusters of lightning strikes and their associated radar echoes. Most lightning 
strike networks have either partial or full weather radar coverage, thus making this 
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validation technique practical for all but a few users. A new technique is presented 
below that uses isolated radar echoes to constrain the error correction and optimization 
procedures to remove the systematic errors. 

Defining the Sample Subset 

A sample subset of lightning strikes is constructed for a short time interval cor- 
responding to low-level radar scans of isolated storms distributed about each DF (refer 
to Figure 14). The interval should account for the propagation of the echo. If the echo 
moves slowly, then the time increment can be increased to enlarge the lightning sample 
size. 


Computing the Bearing Deviations 

Next, the lightning strikes are superimposed on the radar image and the position 
of the reflectivity core is marked. The bearings from each DF to the storm core are 
computed and these are referred to as the "true" bearings. The lightning flashes are 
replotted (Figure 19) and the deviations between the observed and "true" bearings are 
computed. There are 11 flashes associated with this storm in the 10 min interval be- 
tween 2315 and 2325 UTC. 

Flashes 6 and 7 (refer to Table 3) have the two largest semimajor axes of 6.3 km 
and 9.2 km, respectively. Table 5 shows that these are two of the three flashes only 
detected by just two DFs. The low number of flashes seen by DF 2 indicates another 
source of network performance degradation. Poor detection efficiencies and even "blind 
spots" due to site effects (e.g., poor ground conductivity) have been noted by other users 
of such systems . Because one or more sites may not detect a flash, the location must be 
determined from a less favorable geometry (e.g., a flash along the baseline of two DFs or 
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Figure 19. Bearing lines from each of the DFs to the reflectivity core of the 
subject storm. 
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Table 5. DF Bearings and Site Error Corrections 


Flash 

DF 1 

DF 2 

DF 3 

DF 4 

1 

19.9 

5.6 

86.3 

51.4 

2 

18.2 

- 

84.6 

48.9 

3 

19.8 

- 

87.3 

51.3 

4 

18.7 

- 

85.5 

50.7 

5 

19.8 

- 

87.6 

- 

6 

- 

4.1 

- 

51.7 

7 

- 

- 

86.5 

51.2 

8 

19.6 

- 

88.2 

51.1 

9 

18.8 

- 

88.0 

51.9 

10 

19.6 

9.5 

87.2 

51.6 

11 

19.4 

4.6 

88.4 

52.1 

Median 

18.7 

5.1 

87.5 

51.5 

True 

17.1 

350.4 

84.8 

49.0 

Deviation 

-1.6 

-14.7 

- 2.7 

-2.5 
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a flash more distant from one site than another in a more optimal geometry). Over a 
large area, say 300 km from the center of the network, we find that nearly 50% of all 
flashes are only detected by two DFs. 

Because of the small sample sizes, choose the median bearing as the measure of 
central tendency. The median is a more robust estimator in that it is less sensitive to 
outliers or "wild" bearings than is the sample mean. Others have chosen median estima- 
tion for dealing with bearing information for just this reason (e.g., Lenth, 1981). In 
practice, one would try to find a number of storms close to the bearing angles in use 
here (we try to get samples in 6° azimuth bins) to further increase the sample size. The 
median observed bearing and "true" bearing to the storm are also presented in Table 5. 
The nearly 15° deviation at DF 2 is further evidence that the site is poor. Large site er- 
rors such as this in certain directions, however, are not that uncommon (e.g., Schutte et 
al., 1987). 


Computing New Solutions With the Corrected Bearings 

Figure 20 shows a close-up view of the lightning locations before any site error 
corrections are used (shown as dots) and after the initial bearing corrections are imple- 
mented (again described by the 50% error ellipse). The lightning centroid is displaced 
10 km to the southeast of the main echo. Figures 15 and 16 show the lightning locations 
described by the 50% error ellipse after the bearing corrections from Table 5 are imple- 
mented. The lightning centroid is now less than 5 km from the main echo. The impact 
of uncorrected bearing errors on the clustering process and storm identification is dis- 
cussed further in Chapter 4. 
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Figure 20. Location estimates and 50% error ellipses for each flash before 
radar ground truth corrections. 
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Site Error Polynomials 


Once the data base on site error corrections is generated, one would fit a polyno- 
mial to the data having the form 

y = a o + E (sink© + cosk©), k=l,2,3,...,6 (3.1) 

where y is the site error correction to the observed bearing ©. One of the motivations 
for choosing a polynomial fit of this form is that the real-time hardware can implement 
the polynomial (but only to 4th-order) in real-time. Figure 21 shows the site error cor- 
rection curves for the four DFs. The curves were produced by computing site error 
corrections in 6° bins using the technique described by Mach et al. (1986) and Orville 
(1987). A nonlinear regression was performed using the method of Marquardt (1963) to 
compute the polynomial fit to Eqn. (3.1). 

Another way to examine the network performance using the error ellipse is to 
make plots such as Figure 22 where 30 min of lightning data are shown from 2300-2330 
UTC with increasing thresholds of the semimajor axis (SMA). The locations compare 
reasonably well with the radar data even with a semimajor axis of 20 km. 
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Figure 21a. Site error polynomial and residual for DF 1. 
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Figure 21b. Site error polynomial and residual for DF 2. 


56 



Residual (Degrees) Bearing Correction (Degrees) 


DF3 


a a 
p p 


i 


A 

P P 

A A 

P PP P P 


-1 

-2 

-3 

-4 

-5 


* A 

P P 
J 


PP 


P 

A AA A 
P 


P 

A A 


AA 


A 

p P P P 


P 

AAA 

P 

A 

P 


A 
P P 


A A 
P 

A A A 
P P P 


P 

A A 


A A 
P PP 
A 


P P 


P P 
A A 


P P P 
P P 

*0 20 . ... 40. . . . so go.... 100- 120 - • 140- — 160* - -- ItO- - - -200- 220- * - 260- - - - 200- 300* 320- MO - 

e 


A A 


A A A A A 


A A A 

A A A A A 

A 


0 i**-A A — - A -A A A- -A A- — *A«A A A A — -A — -A -A 

U ] AA AAAA A A / A 

! / 

-1 
-2 


.0 20 40 80 *0 100 120.. -.140 100 - • - - 1 * 0 - • . 

e 


.200 2*0 300 320 - * 


Figure 21c. Site error polynomial and residual for DF 3. 
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Figure 2 Id. Site error polynomial and residual for DF 4. 
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Figure 22. Lightning plots during the period 2300-2330 UTC with semimajor 
ellipse axis thresholds: a, infinity; fe, 2 km; £, 10 km; & 20 km. 
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CHAPTER IV. 


CLUSTERING METHODOLOGY 
Storm Identification 

Before generating a lightning time series, a pattern recognition scheme is needed 
to identify the collection of flashes belonging to an individual storm or storm complex. 
For example, Peckham et al. (1984) manually identified storms and storm systems during 
a 3-h sample period using a 2-station lightning network in Florida (Figure 23). Closer 
inspection of Figure 23 shows subareas of greater flash density within the closed con- 
tours representing the boundaries of the 13 storm cells labeled A-M. Given that the 
average thunderstorm lifetime is less than 1 h, these subareas are highly suggestive of 
smaller thunderstorms that existed during a portion of the 3 h observation period. 

In a recent investigation examining the fusion of lightning ground strike infor- 
mation with satellite imagery, Goodman et al. (1988a) subdivided the lightning flashes 
contained within a 4 x 10 s km 2 area into grid cells having a dimension of 0.1° latitude 
by 0.1° longitude. The selected grid cell dimension of approximately 10 km is much 
greater than the random position errors, yet also represents the typical diameter of an 
individual thunderstorm cell. This process was applied during a 1 h period to successive 
S-min sample intervals to permit the manual (human judgment) identification of in- 
dividual and multi-cellular storms using an interactive workstation (Figure 24). The 
closed contours outline the lightning density maxima and can be used to track the move- 
ment, merger, and splitting of clusters. The S-min sample period is commensurate with 
the NEXRAD radar and geostationary weather satellite sampling intervals. In the sec- 
tion below this process is taken a step further and investigate the use of a clustering 
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Figure 24. Cloud-to-ground lightning contour maps (number of discharges per 
cell per 5 min) of storms embedded within a mesoscale convective system in 
North Alabama on 11 June 1986. Time in UTC. (after Goodman et al., 1988). 
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algorithm to assign the individual lightning flashes (the objects) to their respective 
storms (the groups). New members can be assigned iteratively to a new or existing 
storm during each sampling interval and tracked with time. 

Selecting a Clustering Algorithm 

Various clustering methods can be used for such object classifications. Sokal 
and Sneath (1963), Anderberg (1973), Hartigan (1975), and Romesburg (1984), for ex- 
ample, have written entire books on the subject of cluster analysis. At its most basic 
level, clustering is the grouping of similar objects. All clustering algorithms are proce- 
dures for searching through the set of all possible clusterings to find the one that fits the 
data reasonably well (Hartigan, 1975). 

The clustering procedure begins with the choice of some initial partition of the 
data which is then modified so as to obtain a better partition. The basic concept of 
these methods is similar to that of the steepest descent algorithms used for unconstrained 
optimization in non-linear programming (Anderberg, 1973). These algorithms begin 
with an initial point and generate a sequence of moves to find an improved value of the 
objective function until a local optimum is found. The search may involve sorting by 
variables, switching objects between clusters, joining objects together, splitting objects 
apart, adding objects to pre-existing clusters, or specialized searching of a subset of 
clusters (Hartigan, 1975). 

Joining algorithms have been used previously to identify storms from radar 
base-scan reflectivity patterns for the purpose of extrapolating storm motion (Blackmer 
and Duda, 1976; Browning et al., 1982). The basic method follows: 

Step 1 . Identify all gridpoints with local reflectivity maxima above some 
threshold value. 

Step 2 . Compute the spatial separation between each pair of maxima. 
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Step 3 . Combine the nearest pair together to form a new cluster. 


Step 4 . Repeat steps 1-3 until only one cluster remains or a stopping 
condition has been reached (e.g., a threshold on the minimum 
separation distance needed to define a cluster as a storm or storm 
complex). The result can be visualized as a linkage tree or 
dendrogram describing the level (distance) where individual clusters 
join together. 

Step 5 . Repeat steps 1-4 for successive radar images. Compute motion 
vectors from the individual cluster displacements between sampling 
intervals. 


The K-Means Algorithm 


The joining algorithms produce a large number of clusters which are then 
reduced to a manageable number more for the sake of convenience rather than for 
quantitative reasons related to the physical process itself. Chiefly for this reason, and 
because of computational expense (as many as 8000 discharges per hour have been 
detected by the MSFC network), a switching (or transfer) algorithm called the K-means 
or K-nearest neighbor algorithm was selected for grouping the lightning discharges into 
storm clusters. The algorithm begins with an initial partition of the data (referred to as 
seed points) and obtains new partitions until no additional switches in the neighborhood 
of the initial partition improve the classification. Thus, a local rather than a global op- 
timum is sought. The stopping criterion is reached when no movement of an object 
from one cluster to another will reduce the within-cluster total sum of squares. The 
IMSL implementation of algorithm AS 136 developed by Hartigan and Wong (1979) was 
used for this purpose. Appendix C contains the FORTRAN-77 programs developed by 
the author for generating the seeds and calling the IMSL subroutine KMEANS. The 
clustering of the lightning data is implemented as follows: 
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Step 1 . Define an area of interest and subdivide the lightning flashes 
occurring during a sample interval into grid cells having a dimension 
of 10 km x 10 km. Sample intervals of 5-min and 10-min are used 
with the choice a function of either the corresponding radar sampling 
interval (when used for intercomparison or validation experiments) or 
storm lifetime. 

Step 2 . Search the data matrix for isolated lightning density maxima and 
label these as possible seeds. A default threshold of 0.02 discharges 
km' 2 (2 ground discharges per 100 km 2 ) is used almost exclusively 
(see Step 5 below). 

Step 3 . Use the K- means algorithm to obtain optimal partitioning of the 
lightning into storm clusters. 

Step 4 . (For storm tracking purposes). Repeat steps 1-3 for successive 
sample intervals. A storm motion vector can be computed from the 
successive cluster centroid displacements and a time series can be 
constructed from the successive cluster memberships. 

Step 5 . (Hybrid Scheme). Step 2 may generate more seeds than desired 
(or physically realistic) in which case the clusters are not sufficiently 
coherent to track during successive time intervals in Step 4. 
Occasionally too few seeds are generated and Step 3 fails to converge 
to a solution. In the former case with a richness of seeds, the initial 
clusters are joined until the number of remaining lightning clusters 
can be correlated (i.e., tracked) over successive sampling intervals. 

The clusters can be joined further until only a single cluster 
representing an entire complex of storms remains. This final cluster 
can be used to generate and analyze a time series for the entire storm 
system. In the latter case having a deficit of seeds, the seed threshold 
value is lowered to 0.01 discharges km' 2 or a 5 km x 5 km grid 
subarea is searched for additional seeds. 

Step 6 . (Sensor Fusion Scheme). Steps 1-5 address thunderstorm 
identification with the ground discharges alone. An alternative 
seeding method employing radar data has been evaluated. In this case, 
the peak reflectivity tracker or NEXRAD tracker schemes can be 
used to identify the storm cells. These cell centroids then serve as 
the seed points for clustering the ground discharges with the 
K-Means algorithm. 


Demonstration of Method 


Storm systems ranging from individual thunderstorms to large complexes of 
storms are used to evaluate the performance of the clustering procedures. The basic al- 
gorithm is first described for a case of thunderstorm cells embedded within a fast 
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moving pre-frontal squall line that produced widespread severe weather as it crossed the 
Tennessee Valley during a 9 h period in the afternoon of 15 November 1989. The 
widespread precipitation pattern (rain/no rain) and reflectivity peaks (above a threshold 
of 40 dBZ) are depicted at 2000 UTC (1400 CST) in Figure 25. The most active hour of 
cloud-to-ground lightning activity takes place from 2000-2100 UTC. The image is a 
composite constructed from several NWS network (e.g., Nashville, TN; Centerville, AL; 
Jackson, MS) and local warning (e.g., Huntsville, AL; Tupelo, MS) radars in the region. 
An isolated supercell storm ahead of the line (labeled T) was overtaken as the line 
moved northeastward more rapidly than the cell (Figures 26-28). The average system 
motion vector was 25 m s' 1 from 235° (due north is defined as 0°). The merger and 
subsequent interaction of the gust front from the squall line with the supercell led to an 
F4 intensity tornado (estimated wind speeds of 92-116 m s" 1 ) at 2230 UTC that killed 22 
people in Huntsville, AL. The ambient wind in the lower troposphere derived from at- 
mospheric soundings at 1200 UTC (235° at 6.7 m s' 1 at Centerville, AL (CKL) and 230° 
at 8.5 m s' 1 at Nashville, TN (BNA)) is much less than the storm motion vector. In ad- 
dition, the supercell storm had an average storm motion vector (indicated by the large 
dots) of 243° at 18.3 m s' 1 , 8° to the right of the mean wind. 

Define the Region of Interest 

Isolated lightning discharges of both positive and negative polarity occurred both 
ahead of the line from thunderstorm anvils and from the trailing stratiform rain region 
behind the main line of storms. Lightning activity during a 1.5 h period prior to the 
Huntsville, AL tornado is shown in Figure 29. The analysis region of interest is shown 
in Figure 30. The lightning flashes are converted from earth coordinates (latitude, lon- 
gitude) to rectangular coordinates (x,y). For convenience, the geographic coordinate 
(0,0) at the center of the map represents the location of the MIT/Lincoln Laboratory 
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Figure 25. Composite precipitation pattern at 2000 UTC (1400 CST) on 15 
November 1989 observed by the regional network of NWS radars. The radar 
reflectivity is contoured at 18,30.40,45,50,55 dBZ (VIP levels 1-6). 

Small solid circle = Reflectivity cores > 45 DBZ; large solid circle = track 

of the tornadic supercell storm every 15 min. The motion vector of cells within the 

squall line indicated by the wind barb is 25 m s" 1 from 235°. 
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Figure 27. Composite precipitation pattern at 2200 UTC on 15 November 1989. 
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Figure 28. Composite precipitation pattern at 2230 UTC on 15 November 1989. 
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Figure 29. Cloud-to-ground lightning activity from 2105-2236 UTC on 15 
November 1989. ' = negative polarity discharges; + = positive polarity discharges. 
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Figure 30. County and state outlines in the lightning analysis region. Range 
rings are every 50 km. The site of the MIT/Lincoln Laboratory FL2 Doppler 
radar is at the center of the rings. 
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FL2 Doppler radar (refer also to Figure 24). The radar was situated just north of the 
Huntsville airport (HSV) during the COHMEX field campaign to study thunderstorm 
downdrafts, outflows, and gust fronts hazardous to aviation. This transformation is for 
the convenience of some lightning-radar intercomparisons discussed later. Figure 31 
depicts the centroid of all lightning activity in successive 5-min intervals from 2130- 
2235 UTC. The plotted track of the tornadic supercell storm, T, was computed every 
15-min during the interval 2130-2230 UTC from the National Climatic Data Center 16 
mm film-archive of the Nashville, TN (BNA) radar scope. The initial damage from the 
tornado (labeled t) occurred on Redstone Arsenal at 2230 UTC. 

The clustering procedure begins by summing the lightning discharges into 10 km 
grids within an m-row by n-column matrix such as that pictured in Figure 32. The 
evolution of the lightning activity associated with the main line of storms can be fol- 
lowed from the series of contour maps (> 0.02 flashes km" 2 shaded) in Figure 33. A 
contoured lightning map during the 5-min interval 2215-2220 UTC indicates the loca- 
tions of greatest lightning density that might serve as candidate seed points (Figure 34). 
In the following discussion, the gridpoint (20,20) is at the center of the data matrix and 
corresponds to the aforementioned geographic reference point at (0,0). 

Identify the Cluster Seeds 


Step 1 . A copy of the m x n data matrix D is generated and designated 
as the seed matrix S. 

Step 2 . A 3 x 3 point window searches directionally by rows through S 
from the upper left-hand corner (1,1) to the lower right-hand corner 
(m,n) to identify isolated lightning density maxima (Figure 35). 

Any non-zero grid point S(i,j) at the center of the 3x3 window 
having a neighbor greater or equal to itself is also set equal to zero. 
Otherwise, the point is left undisturbed in S as a candidate seed. For 
example, the gridpoint S(15,20)=8 is greater than its neighbors and is 
designated as a seed. When the center of the window reaches the 
gridpoint S(16,20)=6, the value of 6 is less than the neighbor above it 
and is set equal to zero. A single pass through S would produce the 
seeds indicated in Figure 36. However, note that the points 
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Figure 31. Centroid of all lightning within the analysis region in 5-min 
intervals (labeled 1-9-A-D) during the period 2130-2235 UTC. T = track of 
the tornadic storm echo every 15 min beginning at 2130 UTC; 1 = location of 
tornado touchdown and tree damage on Redstone Arsenal. The FL2 radar is 
at the geographic center of the map at the point (0,0). 
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Figure 32. Lightning analysis region with 10 km gridpoints 







2220-2225 


2225-2230 


2230-2235 


Figure 33. Evolution of cloud-to-ground lightning activity in 5-min intervals 
between 2135-2235 UTC. Contours are approximately 0.01 discharges km' 2 
and 0.02 discharges km* 2 (shaded). 
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Figure 34. Contoured lightning density map during the period 2215-2220 UTC. 
Contour interval is every 0.01 discharges km 2 beginning with the value 0.02 
discharges km' 2 . 
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Figure 35. Diagram of the 3x3 point search window moving through the 
m-row x n-column data matrix D. Lightning shown for the period 
2215-2220 UTC at each 10 km gridpoint. 
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S(13,22)=2 and S(17,20)=2 are not isolated maxima in D, but result 
from the left to right search direction of the 3x3 window. Such 
spurious seeds are removed by step 3. 

Step 3 . The original data matrix D is scanned locally for a neighbor that 
might be larger than the candidate seed itself, but had been set to 
zero in S in the prior step. If a neighbor is not a seed, but is still 
greater than the candidate point in question, then the candidate point 
is not an isolated maxima and is set to zero. For example, let the 
window continue moving through S until the center of the window is 
at S(17,20)=2. S(17,20) is now a candidate because S(16,20)=0 
resulted from step 2. However, S(17,20) < D(16,20) so S( 17,20) 
should also be set to zero. Thus, one pass produces the final S 
matrix of exactly K seeds (Figure 37). 

Step 4. Define a set of criteria for establishing the number of valid 
seeds. In this example, there are 17 seed points having at least 0.01 
flashes km' 2 (or 1 flash within a 10 km x 10 km grid). The possible 
seeds consist of 7 seeds produced by single lightning discharges of 
positive polarity, 5 produced by single discharges of negative 
polarity, and 5 seeds with a density greater than 0.02 flashes km A 
sensitivity study of possible criteria and their justification are 
discussed below. 


Determine the Cluster Memberships 


The initial partition of K seeds is critical since the K-means algorithm must op- 
timally assign the lightning flashes to one of exactly K clusters. A different initial par- 
tition (determined by the number and location of the seeds) might produce a different 
final partition and storm motion vector based on the tracking of the seeds. Furthermore, 
any change in cluster membership also alters its time series. 

Figure 38 shows the the cluster assignment for each flash occurring between 
2215-2220 UTC and the location of the 5 seeds (labeled a-e) identified in step 4 above 
with a density of 0.02 flashes km' 2 or greater. Seed densities of 0.01 flashes km' 2 are 
excepted and the location of the tornadic supercell storm (labeled T) as observed by the 
BNA radar is indicated. The Tennessee-Alabama border is approximately 45 km north 
of the reference point. 
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Diagram of the seed matrix S after one pass through the data matrix, 
window is centerd at the point S(i,j) = (17,20). 



Column 


Figure 37. Diagram of the final seed matrix S. 
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Figure 38. Cluster assignments for each flash during the period 2215-2220 UTC. 
Circled lower case letters fa-e) = cluster seeds with values greater than 0.02 
discharges km' 2 ; T = location of tornadic storm echo at 2215 UTC. 





Clusters B,C, and D have the smallest displacements between their seed (first 
guess) and final centroid. Clusters A and E are more spread out with the four 
northwesternmost flashes assigned to cluster A being of positive polarity. A comparison 
with the lightning density contour maps (Figures 33; 34) and the regional radar image at 
2215 UTC (Figure 27) indicates good agreement between the 5 lightning centroids and 
the largest thunderstorm echoes. The positive polarity flashes in the northwest quadrant 
of cluster A are seen to be located in the trailing stratiform rain behind the main line of 
thunderstorms. 

Next, let the 5 single discharges of negative polarity serve as additional seeds 
(f-j) and consider the impact on the preceding example (Figure 39). The distribution of 
the 10 seeds splits apart cluster E into E, I, and J. Cluster A is split into clusters A and 
F, where F is comprised entirely of the scattered positive polarity flashes in the trailing 
stratiform rain. Cluster B is subdivided into clusters B and G and lastly, cluster H is 
split off from cluster C. Comparison with the 2215 UTC radar image indicates that 
cluster H is probably best represented by the solitary echo southeast of the main line of 
thunderstorms. However, the remaining subdivisions appear less realistic. 

The effect of changing cluster memberships is presented in Table 6. The original 
cluster membership changes (and decreases) when the number of seeds increases from 5 
to 10 seeds. Cluster D maintains the same total membership, but exchanges members 
with clusters C and E. In addition, the within group sum of squares (WSS) decreases as 
the number of seeds is increased. 

Due to the natural spatial variability of the lightning strikes, a single flash within 
a 100 km 2 grid is just as likely to have occurred in any of the bordering grids. Indeed, 
ground-based radar and high altitude airplane lightning measurements of large weather 
systems such as this have shown lightning flashes occasionally propagating over 100 km 
horizontally before striking the earth (Ligda, 1956; Goodman et al., 1988c). 
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Figure 39. Cluster assignments for each flash during the period 2215-2220 UTC 
Pireled lower case letters (a-e) = cluster seeds with values greater than 0.02 
flicharges ^~ z ~ Ckskd jgw£L case letters (f-i) = cluster seeds each 
represented by a single negative polarity discharge. 




Table 6. Cluster Characteristics at 2215-2220 UTC 



84 


Tracking Seeds and Clusters 


The eastward propagation of the lightning activity can be followed over the three 
subsequent 5-min periods 2220-2225 UTC (Figure 40), 2225-2230 UTC (Figure 41), and 
2230-2235 UTC (Figure 42). In each case a threshold seed density of 0.02 flashes km* 2 
is employed to identify the clusters. The seed and cluster letter identifiers are newly as- 
signed each 5-min sample period in the order the seeds are identified. Thus, the assign- 
ment follows the search direction from upper left to lower right within the seed matrix 
S. 

Cluster A beginning at 2215 UTC can be tracked from the continuity between 
5-min observation periods from 2215-2230 UTC as A- A- A- A. Since no local seed is 
identified at 2230 UTC, the singular A event near the coordinate (70,120) in Figure 42 
is misclassified as a member of the cluster to its south, rather than to the decaying 
northern storm. Cluster B beginning at 2215 UTC can be tracked as B-(B+C)-B-A. The 
tornadic storm cluster can be tracked as C-(D+E)-C-C. The two clusters initially iden- 
tified at 2215 UTC as B and C are each split apart into two distinct groups at 2220 UTC 
before merging again. The additional lightning clusters at 2220 UTC makes subsequent 
tracking of all but the largest storms difficult. However, such storms also produce a 
greater number of discharges and hence pose the greatest threat. Isolated storms, large 
or small, tend to maintain their identity. Storms that merge and split apart cause iden- 
tification problems for the radar echo tracking techniques as well. 

An approach that might be considered for reducing the number of clusters to 
only those that are coherent (i.e., trackable) from sample to sample is to make use of a 
hybrid scheme using multiple methods. One possible hybrid approach has been applied 
to this problem. The K-means algorithm is run initially as before, but then a variation 
of the joining algorithm is used to reduce the number of clusters to the same number 
identified in the prior sampling interval. The final number of clusters desired can be 
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Figure 40. Cluster assignments for each flash during the period 2220-2225 UTC. 
Circled lower case letters (a-i) = cluster seeds with values greater than 0.02 
discharges km' 2 . 
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Figure 42. Cluster assignments for each flash during the period 2230-2235 UTC. 
Circled lower case letters (a-O * cluster seeds with values greater than 0.02 
discharges km’ 2 . 
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based on having a manageable number of entities to track, correspondence with radar 
measurements, or perhaps based on some criteria such as desired storm size (by setting a 
minimum areal extent threshold). The result is shown in Figure 43 for the period 
2220-2225 UTC. The hybrid implementation involves finding the mean (x,y) coordinate 
of the two nearest clusters to be paired until the total number of seeds is reduced from 
the original set of 9 seeds (refer to Figure 40) to a revised set of 5 seeds. Table 7 lists 
the results of changing the the number of seeds. When clusters D and E are joined, the 
seed location of cluster E is used instead of a mean location for the pair because of the 
greater seed density and greater number of members assigned to E in the first applica- 
tion of the K-means algorithm. 

Alternatively, the grid cell dimension could be increased to 20 km x 20 km or 
greater to allow greater seed densities and fewer seeds, but this approach will reduce the 
accuracy of the seed locations and the ability to resolve small storms as individual en- 
tities. We have found that it is better to generate more seeds and iteratively join the 
clusters than to generate too few seeds and be unable to resolve new cells that may 
produce low flash rates or identify decaying storms which produce scattered discharges. 
The Huntsville, AL tornadic storm, for example, produces a flash density of 1 discharge 
per 100 km 2 earlier in its life-cycle. Other alternative approaches using only lightning 
data, considered to be beyond the scope of this study, would be 1) to set a threshold on 
the seed density or total cluster membership at some level greater than 0.02 flashes km 
that might be considered to be physically meaningful, or 2) define a maximum search 
radius, say 20 km, which would be the maximum distance a flash could be from a seed 
point for consideration as a cluster member. 

Figure 44 shows 30-min seed tracks (representing the gridpoint with the largest 
local flash density) from 2200-2230 UTC for the clusters identified at 2215 UTC as 
A-E. The lower case letters identify only the 5-min periods when the seed value is 
greater than the 0.02 flashes km' 2 threshold. The upper case letters indicate the location 
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Figure 43. Hybrid clustering algorithm assignments during the period 2220-2225 UTC. 
Circled lower case l etters (a-e) = cluster seeds with values greater than 0.02 
discharges km' 3 . 
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Table 7. Hybrid Clustering Algorithm at 2220-2225 UTC 
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Figure 44. Seed tracks during the period 2200-2230 UTC. Lower case letters 
“ seed position for each 5-min interval that the seed value exceeds a threshold 
of 0.02 discharges km' 2 . Upper case le tters tA-El = seed seed location at 
2215 UTC; T = Position of tornadic storm echo at 2200, 2215 and 2230 UTC. 
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of the seeds at 2215 UTC and the location of the tornadic storm echo at 2200 UTC, 
2215 UTC, and 2230 UTC. All clusters are generally moving from west to east, with 
lightning cluster C within 1-gridpoint of the storm echo. However, only seeds B and D 
are greater than the threshold value during each 5-min interval. The track variations 
suggest a moving or weighted average of past tracks be used instead of updating the 
track with each observation. Such techniques are also employed in the radar echo track- 
ers described earlier. 


Synergism W ith Radar 


A primary objective of using a clustering algorithm for pattern recognition is to 
allow the identification and tracking of storms and their changing membership using the 
lightning data alone. However, the previous intercomparisons between the lightning 
clusters and the radar echoes suggest that a fusion of the data sources would be useful. 
For example, the radar echo trackers could provide seeds for clustering the lightning. 
Consider the tornadic storm at 2135 UTC when the maximum flash density is only 1 
flash per gridpoint. Using a seed threshold of 0.01 flashes km 2 produces 21 total 
clusters and permits a separate cluster of two flashes, labeled T, for the tornadic storm 
(Figure 45). Using a seed threshold of 0.02 flashes km causes the two lightning flashes 
to be assigned to cluster C. 

On the other hand, the size and intensity criteria used in radar echo tracking 
could be given adaptive thresholds to capture small, electrically active storms that do not 
meet the default threshold criteria. In a study of the effect of radar echo size and 
asymmetry on the identification of small, electrically active microbursts occurring in 
Alabama and Florida, Buechler and Goodman (1990) find that the operational NEXRAD 
storm identification algorithms have a probability of detection less than 0.5. The algo- 
rithms are designed to detect large severe storms. Small storms (less than 5 km in 
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Figure 45. Cluster assignments during the period 2135-2140 UTC. A seed 
threshold of 0.01 discharges km' 2 produces 21 clusters (A-T). X ” tornadic 
storm discharges. 
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diameter), changing echo shapes, and mergers sometimes cause the radar echo trackers to 
lose (fail to identify) a storm from one 5-min scan period to the next. The continued 
occurrence of lightning could thus be used to alert the echo tracker to adjust threshold. 

A comparison of the different lightning clusters produced from radar echo seeds 
follows. The NEXRAD storm identification and tracking algorithm (Bjerkaas and For- 
syth, 1980) is implemented with S-band radar data collected on 13 July 1986 at 2328 
UTC by the CP2 radar. The tracking algorithm yields two echo centroids, labeled A and 
B, for this multi-cellular complex of storms (Figure 46). A smaller length threshold (5 
km) and intensity threshold (30 dBZ) could produce more than two distinct storms. A 
peak reflectivity tracker (e.g.. Crane, 1979; Rosenfeld, 1987) would produce at least four 
echo centroids, indicated by the shaded area representing reflectivity values in excess of 
55 dBZ. The 3-dimensional structure of these storms is portrayed in Figure 47. The 
lightning seeding algorithm (with its 10 km grid) would produce only a single seed lo- 
cated between the points (36,115) and (46,125) depending on the placement of the grid. 
In this isolated storm example, a joining-type clustering algorithm could be used instead 
to assign the 28 lightning flashes to possibly four distinct lightning clusters. The lightn- 
ing locations in Figures 48 and 49 suggest four storm clusters offset ahead of the peak 
reflectivity maxima. The offset is partially physical in that the lightning strikes 
generally occur outside of the storm reflectivity core, but there also appears to be a 
directional bias to the southeast of the echo centroids, likely due to unrecovered sys- 
tematic errors at one or more of the antennas. This initial examination of the perfor- 
mance and limitations of lightning clustering algorithms demonstrates their usefulness in 
identifying and tracking thunderstorms. In addition, it appears that more than one 
method and stopping criteria is best suited to the various space and time scales examined 
here. 
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Figure 46. Plan-view of a small multicellular storm complex observed by the 
CP2 radar on 13 July 1986 at 2328 UTC. Circled upper case letters CA-BI = 
two storm centroids identified by the NEXRAD storm identification and 
tracking algorithm; solid contour - reflectivity > 40 dBZ; shaded region - 
reflectivity > 55 dBZ. Distance units in kilometers from CP2. 
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Figure 47. Three-dimensional structure of the storm complex in Figure 46. 
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Figure 48. Cluster assignments using two centroids identified by NEXRAD 
algorithm. Circled l ower case letters (a.b) = radar echo cluster seeds; 

Upper case letters (A.B1 * cluster assignments. 
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CHAPTER V. 


THUNDERSTORM LIFE-CYCLES AND EXTRAPOLATION FORECASTING 
Testing and Evaluation of the Logistic Growth Model 

The storm identification and tracking process using clustering methods sets the 
stage for assessing the validity of the logistic growth model and predictability of lightn- 
ing activity during the storm life-cycle. Table 8 gives the model parameters for a 
selected number of cases. A total of 14 storm systems have been examined ranging in 
size from small, short-lived air mass storms (< 1 h duration) to large, long-lived mesos- 
cale convective systems. The total (intracloud and cloud-to-ground) lightning life-cycle 
has been computed for 3 of the short-lived storms. These latter cases are included to 
demonstrate the applicability of the logistic model to storms for which the total lightning 
rates portray the growth and decay process, yet the storm may produce too few discrete 
cloud-to-ground flashes to establish a trend (e.g., the 20 July case). The cloud-to- 
ground lightning data were summed over 1, 5, 10, 30 and 60 min sampling intervals 
depending on the availability of data or the duration of the storm life-cycle. The SAS 
(1985) non-linear regression procedure NLIN using the Gauss-Newton method was used 
to estimate the parameters /9 and k from each set of observations and limiting values of a 
corresponding to the total number of lightning flashes produced by each storm event. 

Isolated and Multicellular Storms 

Figure 50 shows the evolution of the cloud-to-ground lightning activity as- 
sociated with 3 distinct storms observed in Southern Tennessee on 17 July 1986. The 
storms were identified from the cluster analysis during the interval 1605-1845 UTC 
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Figure 50. Cloud-to-ground lightning clusters during 5-min intervals on 17 
July 1986. Distance units in kilometers from FL2 radar. (Time in UTC.) 
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(1105-1345 CDT). The maps were produced every 10 min from 1605-1735 UTC, a 
period which encompasses the entire 90 min life-cycle of cell A. The corresponding 
radar echo time history is summarized in Figure 51 from 30 min observations made with 
the Nashville, TN radar. The center coordinate of each radar map is approximately 6 
km east and 9 km north of the FL2 radar (used as the center of the lightning maps). 

The lightning and rainflux time series for storms A and C are shown in Figure 
52. Cell A undergoes a secondary surge in lightning activity after an initial peak. The 
rainflux time history also shows two peaks, each lagging the lightning maxima by 10 
min. Cell C shows a single maxima in both lightning activity and rainflux. 

The logistic model regression (P), observed flash rate (A), and residual error (A- 
P), for storms A and C are plotted in Figures 53 and 54. The model fit is excellent in 
both cases with correlation coefficient r > 0.998. Cell A produced 102 ground discharges 
(the observed a) and Cell C produced 453 ground discharges during its 120 min life- 
cycle. The steepness of the curves along the time axis are described by the parameter k. 
For storm lifetimes of 1-2 h (sampled at 10 min intervals), the median value of k (Table 
8) is approximately 0.8. 

Total Lightning Activity for Airmass Storms 

The value of k increases to 0.9 for storms with lifetimes less than about 1 h. 
Figure 55 shows the logistic model results for the total lightning produced by the 20 July 
airmass storm described in Chapter 2. Due to the short-lived 13 min duration of the 
lightning activity a more frequent sampling period of 1 min was chosen. Again, the 
residual error is only a few percent. The other two storms occurred in Florida in 1977 
and 1978 (Figure 56), and have been studied previously by Piepgrass et al. (1982) and 
Krehbiel (1986). The 11 July 1978 single peak storm is similar to the 20 July storm, but 
the secondary peak of the 8 August storm is likened more to Storm A on 17 July 1986. 
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Figure 51. Base scan radar echoes observed by the Nashville radar upper left . 
1600 UTC; upper right . 1630 UTC; lower left . 1700 UTC; lower right . 

1730 UTC. Contour interval every 10 dBZ beginning at 10 dBZ. 
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Figure 53. Logistic model regression and residual error for 17 July 1986 Storm A. 
Top: For each 10 min observation period, & - the observed flash rate; 

P = model prediction. Bottom: Residual error - (A - P). 
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Figure 54. Logistic model regression and residual error for 17 July 1986 Storm C. 
Too: For each 10 min observation period, A = the observed flash rate; 

P = model prediction. Bottom: Residual error = (A - P). 
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Figure 55. Logistic model regression and residual error for 20 July 1986 airmass 
storm. Tod: For each 10 min observation period, A = the observed flash rate; 

P * model prediction. Bottom: Residual error = (A - P). 


108 



Total Flash Rate (per 5 min) 



Figure 56. Total lightning time series (discharges per 5 min) for two thunderstorms 
observed near Cape Canaveral, Florida. 
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In the latter case the total system life-cycle is nearly symmetric within the envelope en- 
compassing the storm, but the large individual peaks indicate 2 or 3 embedded cells or 
resurgent growth and decay. Few ground discharges were observed in all three cases. 

Long-lived Storm Complexes 

Figure 57 depicts the life-cycle of the 15 November 1989 mesoscale weather sys- 
tem through much of its life-cycle. The time series covers the period 1600 UTC 15 
November to 0100 UTC 16 November in 10 min increments. Despite the multiple peaks 
superimposed upon the curve (due primarily to cell mergers), the envelope of the lightn- 
ing activity exhibits symmetry about the maxima of 360 flashes which occurs at time 
period 27 (2030-2040 UTC). There are 19 time steps between 100 flashes and the maxi- 
mum during the growth phase and 17 time steps from the maximum to 100 flashes 
during the decay phase. This 2 time step difference is just 20 min over a 6 h period. 

The logistic growth curve in Figure 58 again fits the data well, thus reinforcing 
the concept of symmetric growth and decay at yet larger space and time scales. The 
long-lived (>2 h) storm systems have k values less than about 0.4, with k inversely 
proportional to storm lifetime (Table 8). Using the symmetry that characterizes logistic 
growth, we can estimate that the duration of the decay phase will approximately equal 
the duration of the growth phase. This estimate also defines the valid extrapolation 
period for yes/no (presence/absence of lightning activity) forecasts. Forecasts of the 
duration of these large storms are valuable because their long lifetimes and severe 
weather production makes them the most disruptive and hazardous weather systems 
during the Spring and Summer. 

Figure 59 shows ten individual lightning ground strike time histories and a com- 
posite life-cycle for convective storm complexes observed in the Oklahoma/Kansas 
region of the Southern Great Plains with the NSSL ground strike network (Goodman and 
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Figure 58. Logistic model regression and residual error for 15 November 1989 
mesoscale weather system. Too : For each 10 min observation period, A = the 
observed flash rate; £ - model prediction. Bottom: Residual error - (A - P). 
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Figure 59. Cloud-to-ground lightning discharge histograms and composite lighting 
life-cycle for mesoscale convective complexes (MCCs). The four life-cycle phases 
are identified as first storms (F), initiation (1), cold cloud shield maximum extent (M), 
and termination (T). The composite is plotted with respect to the time and 
magnitude (± one standard deviation) of the average peak flash rate. 

(After Goodman and MacGorman, 1986). 
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MacGorman, 1986). The life-cycle is identified from the infrared cloud temperature 
criteria developed by Maddox (1980). The dimension of the cold cloud top observed by 
satellite is used to define four phases of the storm system life-cycle. These phases iden- 
tify the formation of the first storms (F), the initial time at which the size criteria 
threshold is reached (I), maximum extent of the cloud shield (M), and the final time at 
which the size threshold is exceeded (T). 

In Figure 60 the composite life-cycle has been normalized to the maximum 
hourly flash rate and to the time at which it occurred. The flash rates increase and 
decrease exponentially with time. The exponential relations best fitting the data are also 
shown where N is the fraction of discharges in a given hour occurring at a time, t, rela- 
tive to the magnitude and occurrence of the peak. The logistic model fit to the com- 
posite is given in Figure 61. 


Interpretation of Results 

In general, the longer a convective storm takes to reach its maximum intensity, 
the longer it takes to decay. A storm cell that grows rapidly often decays rapidly. The 
logistic model provides a good fit to the observed data. In each case the model residual 
error is only a few percent for individual time series of cloud-to-ground lightning and 
total lightning activity. The rate constant k is greatest for the shorter-lived storms (refer 
to Table 8). The limiting value, a, is a function of storm duration and the environmen- 
tal factors that serve to sustain strong convection (e.g., atmospheric instability, a mois- 
ture source, and a triggering mechanism such as an upper level jet streak, front, storm 
merger, or outflow boundary). The limiting value a and residual error will now be con- 
sidered in greater detail. 
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MCC LIGHTNING COMPOSITE 



Figure 60. Average hourly cloud-to-ground discharge rates are normalized to the 
average peak flash rate of MCCs and are shown relative to the time of occurrence 
of the peak. N is the percentage of the peak rate at a given hour and R is the 
correlation coefficient for each of the exponential curves. Open circles, denote 
the time and magnitude of the lightning rates for each of the MCC life-cycle 
phases F, I, M, T (After Goodman and MacGorman, 1986). 
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Figure 61. Logistic model regression and residual error for MCC composite 
life-cycle. Too: For each hourly observation period, A = the observed flash rate; 
E « model prediction. Bottom : Residual error = (A - P). 
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Role of the Environment in Limiting a 


One frequently used index to describe the instability of the storm environment is 
the lifted index. The lifted index (°C) is found by lifting a positively buoyant saturated 
parcel of air, having the thermodynamic characteristics of the lowest 100 mb layer of 
the atmosphere, from its level of free convection (where the buoyancy is first positive) 
to 500 mb. The temperature excess of the parcel,T y , to that of the environment, T, at 
500 mb is the lifted index in °C. Figure 62 shows the relationship between maximum 
hourly flash rate during June 1986 within the 8 state area shown in Figure 29 and the 
1200 UTC lifted index computed from soundings made at Redstone Arsenal, AL. Total 
flash rates are seen to increase non-linearly with decreasing atmospheric stability. 

A related parameter that describes the environment and offers additional physical 
insight is the convective available potential energy (CAPE). CAPE is proportional to the 
square of the maximum parcel updraft speed, w 2 , and thus represents the increase in 
kinetic energy of a parcel associated with its vertical acceleration (Weisman and Klemp, 
1982). CAPE is also commonly referred to as the available buoyant energy, i.e., the 
kinetic energy due to buoyancy. Figures 63 and 64 show, respectively, the relationship 
between 1) total cloud-to-ground lightning and CAPE; and 2) maximum rain rate and 
maximum parcel energy (CAPE as defined by Zawadzki et al., 1981). Both cloud-to- 
ground lightning (r=0.62) and rainfall (r=0.79) are positively correlated with the buoyant 
energy in the environment. It should be noted that the predicted value of w ignores the 
effects of precipitation (or mass) loading, vertical pressure gradient perturbations, and 
mixing which would lower w estimates by as much as 50%. Yet, it is suspected that the 
total lightning rates would produce an even greater correlation. Such measurements will 
not be possible for large storm systems until total lightning observations are available 
from space in the late 1990s with NASAs lightning mapping sensors (see Chapter 7). 
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12 UTC Lifted Index (°C) 


Figure 62. Regression plot of maximum hourly flash rate (Y) as a function of the 
lifted index (X) computed from the Redstone Arsenal 1200 UTC soundings taken 
during June 1986. The 95 percent confidence limits are indicated by the dotted lines. 
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Figure 63. Regression plot of total cloud -to -ground flashes observed each day (y) 
is a function of Convective Available Potential Energy (CAPE) (x) computed from 
he Redstone Arsenal 1200 UTC soundings taken during June and July 1986. 

The 95 percent confidence limits are indicated by the dotted lines. 
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Figure 64. Maximum rain rate as a function of maximum parcel energy for 67 
storm days near Ottowa, Canada during the summer of 1969-1970. 

(After Zawadzki et al., 1981). 
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An Examination of the Residual Errors 


Although we claim the logistic model describes the thunderstorm life-cycle 
reasonably well, there may be other non-linear models that describe the life-cycle just as 
well, if not better. The logistic model offers a high degree of correlation with the ob- 
served data and errors are small, but we note that the residuals are consistently above or 
below the fitted values for short periods. Such patterns suggest positive autocorrelation 
in the error terms. Such correlation patterns further suggest improved models can be 
constructed by adding one or more independent variables to the present model. We can 
evaluate the error terms quantitatively with the Durbin-Watson test statistic (Neter et al., 
1983). The test statistic, D, is defined as 


D = 



(5.1) 


where n is the number of sample periods, e t is the residual error from the least squares 
regression, E e t 2 is the residual sum of squares, and the term (e t - e t l ) is the difference 
in the residuals at two successive times. The usual hypothesis test alternatives are 


H : p = 0 

O r 

Hj : p > 0, 


(5.2) 


where p is the autocorrelation parameter. The decision rule is as follows: 


If D > dy, conclude H o , 
If D < d L , conclude H r 


(5.3) 
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where d u and d L are the upper and lower bounds obtained by Watson and Durbin, such 
that a value of D outside these bounds leads to a decision and any d L < D < d w gives an 
inconclusive result. 

Estimating the two parameters f) and k (a is proscribed), results in d L -1.34 and 
dy-1.42 at a level of significance a o -0.01 (Table A-6 in Neter et al., 1983). For ex- 
ample, the 15 November 1989 storm system has n-53 and a value of D-0.06 « d L . 
Thus, one concludes the error terms are positively correlated. For the MCC composite 
n-18, d L «0.90, and dy-1.12, D-0.84 < d L , and again conclude H r Small values of D 
usually lead to the conclusion that p> 0 because successive error terms tend to be of the 
same magnitude. The residual errors tend be largest when cells merge or cells 

redevelop/intensify, both factors which lead to a short term increase in the observed 
flash rates. The 15 November 1989 case shows a number of such peaks associated with 
cell mergers. Much of this variation is reduced in the MCC composite by averaging the 
ten cases. These smaller scale interactions can be included in the basic logistic model by 
the inclusion of additional (e.g., exponential) terms. 

Applications to Extrapolation Forecasting 
Thunderstorm Duration 

The relationship between the environmental instability and energy can be used as 
a "first guess" of a. The storm system duration will increase and the system growth rate 
will decrease as a increases Thus, a can serve as an index to provide some insight into 
the potential life-cycle of the storms to the forecaster. Since the logistic curve is sym- 
metric about its inflection point, knowledge of when a storm or storm complex has 
begun the decay phase of its life-cycle (as indicated by decreasing lightning discharge 
rates) can be used to estimate the valid extrapolation range. One might also predict that 
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lightning activity will decrease exponentially in the same time it took to reach its peak. 
This type of information can be used in the context of yes/no (lightning/no lightning) 
forecasts. The existing extrapolation forecasting procedures use an arbitrary forecast 
period determined by some ad-hoc method. Using an appropriate storm or system mo- 
tion vector, one now has a physical basis and empirical evidence upon which to make 
the determination of longer or shorter extrapolations. 

Thunderstorm Existence 

Consider the yes/no forecast of lightning activity for the 15 November 1989 
storm system using the lightning grids. A total of 1269 cloud-to-ground flashes oc- 
curred in association with the line of storms that came through the Tennessee Valley 
during the 1 h interval 2130-2235 UTC within the analysis sub-region shown in Figure 
33. The time series (in 5-min intervals) over this 1 h period is shown in Figure 65. 

A full intensity extrapolation forecast would simply move the lightning in one 
map with the storm motion vector At (5-min) steps forward. If the storm motion vector 
placed the system and individual embedded storms in the correct location t+At step 
ahead, the existence question would be answered correctly. 

Since the total cloud-to-ground lightning for the system is changing (growth and 
decay), the forecast error will depend on the number of steps ahead to forecast. For ex- 
ample, the system motion vector of the lightning from 2130-2135 UTC to the next 
5-min interval is 4.8 km from the southwest at 246° (due north is 0°). An extrapolation 
of the pattern t+At steps ahead is simply 

x * x o + nAt(Ax) (5.4) 

y = y o + nAt(Ay) (5-5) 
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cloud " to “8 round lightning time series (per 5 min) of the 15 November 
1989 storm system during the period 2130-2235 UTC. 
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where n is the number of sample periods ahead, and Ax and Ay are the x and y com- 
ponents of the storm motion vector. The displacement of each lightning discharge is 
computed first and then a grid is generated at the desired time. Consider the lightning 
activity at 2130-2135 UTC which is extrapolated in space 11 periods ahead (Figure 66). 
Compare this forecast grid against the observed grid at 2225-2230 UTC (Figure 67). A 
qualitative evaluation of the forecast can be made using the evaluation criteria developed 
by Donaldson et al. (1975) and successfully used by Browning et al. (1982) and many 
others for evaluating radar echo extrapolations. The scoring criteria for lightning are 
given in Table 9. The evaluation criteria of interest are the Threat Score or Critical 
Success Index (CSI), probability of detection (POD), and false alarm rate (FAR) defined 


CSI = 

A 

(5.6) 

(A+B+C) 

POD = 

A 

(5.7) 

(A+B) 


FAR = 

C 

(5.8) 

(A+C) 


A is defined as the number of gridpoints predicted to have a flash density >0.02 
km' 2 at the valid forecast time (t+At) and an observed density of 0.01 km' 2 or greater 
within a local neighborhood of the gridpoint. A 3x3 point verification grid centered at 
the gridpoint of interest is used for this purpose. B is defined as the number of grid- 
points where there is no flash density >0.02 km' 2 predicted at the valid forecast time, 
yet there is a gridpoint observed in its neighborhood with a density >0.02 km' 2 . Lastly, 
C is defined as the number of gridpoints predicted to have a flash density >0.02 km , 
yet no lightning is observed in the local neighborhood. The validation focuses on hitting 
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Figure 66. Eleven-period ahead extrapolation (valid 2225-2230 UTC) of lightning 
activity observed during the period 2130-2135 UTC. 



Figure 67. Observed lightning activity during the period 2225-2230 UTC. 
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Table 9. Objective Forecast Scoring Criteria 


Observed 

Forecast 

Lightning 
No Lightning 

Lightning 

No Lightning 

A 

B 

C 

D 










or missing the major storm centers and the validation area is similar to that of Browning 
et al. (1982) who used 20 km x 20 km grids to validate their rain/no rain extrapolation 
forecasts. 

The 11 period ahead forecast produced a CSI-0.46, POD-0.71, and FAR-0.43. 
The residual error e t =3 for forecasting the number of seeds or clusters. A total of 8 
seeds were identified at 2130-2135 UTC (forecast to still be extant at 2225-2230 UTC 
since the full system was moved forward in time) and only 5 seeds are observed (refer to 
Chapter 4). The loss of storms most severely reduces the CSI and FAR since we assume 
the storm clusters will still be present nearly an hour later. Due to the birth and death 
process, the CSI and FAR tend to improve with decreasing lag. We note that the lightn- 
ing cluster identified as I (Figure 45) has long since dissipated (2140 UTC), yet it is the 
forecast valid at 2225-2230 UTC. 

Lastly, consider the 1 -period ahead extrapolation of the lightning activity from 
(2225-2230 UTC) to (2230-2235 UTC) (Figure 68) and compare this grid with the ob- 
served lightning grid at (2230-2235 UTC) (Figure 69). We find the CSI-0.67, POD=0.71, 
the FAR=0.08, and the number of clusters increases from 5 to 6. In this example the 
FAR is greatly reduced with a lesser improvement in the CSI. Although we generally 
expect better forecasts at shorter lags, a number of varied weather scenarios need to be 
studied to understand ways to improve the forecast scores. 
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Figure 68. One-period ahead extrapolation (valid 2230-2235 UTC) of lightning 
activity observed during the period 2225-2230 UTC. 
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Figure 69. Observed lightning activity during the period 2230-2235 UTC 
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CHAPTER VL 


SUMMARY AND CONCLUSIONS 

The goal of this research was to develop a mechanistic model that can be used 
for extrapolative forecasting of lightning and thunderstorm activity. The original 
research reported herein has resulted in the development of a 3 parameter logistic 
growth model to explain the exponential growth and decay of lightning activity accom- 
panying thunderstorms, at different space and time scales extending over four orders of 
magnitude. The logistic model depicts a process where the growth rate is proportional to 
the size reached, the remaining size of the population, and is a function of time. The 
growth rate constant depends on the size of the storm while the limiting value of the to- 
tal lightning activity and lifetime is related to the available energy in the environment. 

Short-lived storms may not produce sufficient cloud-to-ground lightning dis- 
charges to exhibit a continuous growth and decay life-cycle. Instead, the total 
(intracloud and ground) discharge rate best describes the thunderstorm life-cycle. In 
two severe weather cases, a short-lived microburst-producing storm and a long-lived 
tornadic supercell storm, peak total lightning rates were greater than 22 flashes min' 1 
and 40 flashes min' 1 , respectively. Yet cloud-to-ground flashes occurred at rates of 
only 1-5 min' 1 . The logistic growth model using cloud-to-ground lightning data is more 
appropriate for describing the evolution of storm complexes, where the longer sampling 
period is less affected by sudden surges in growth. The short-lived cells produce too 
few samples (ground discharges) to adequately describe the system and detect significant 
trends. 
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The model residual errors, though small, have been shown to exhibit positive 
serial correlation. This is a result of the merger and reinvigoration of storms giving rise 
to multiple local maxima in the lightning time series during the growth and decay of the 
storm system. Suggested techniques for resolving or deconvolving these multiple peaks 
to gain further insight into their behavior (and implications for nowcasting) might in- 
clude other non-linear models with additional parameters, stochastic time series models 
(Box and Jenkins, 1976; Abraham and Ledolter, 1983), and simple first-order difference 
equation formulations of the logistic model that are used to describe the chaotic behavior 
of non-linear dynamical systems. The predator-prey relationships in animal and biologi- 
cal populations have been modeled in this way (May, 1974; May 1976). 

Many of the physical (non-instrumental) sources of forecast error examined have 
the same root causes that affect the radar echo extrapolation systems, i.e., incorrect ex- 
trapolation vectors, new storms (births), renewed growth (due to cell mergers or a 
change in the storm environment, storm decay or total dissipation (deaths). The most 
promising improvements to extrapolation forecasts should come from the use of ap- 
propriate non-linear models to understand the birth and death process. It appears ap- 
propriate to maintain a history on the life-cycle of the system when two or more storms 
merge, and not, as current methods do, begin a new description from the time of the 
merger. This will lead to incomplete life-cycle descriptions and delete the information 
needed to predict the decay of storm complexes. 

A novel constrained optimization approach was developed to remove systematic 
errors from the cloud-to-ground lightning data base. An optimization algorithm con- 
strained by the observed position of isolated radar echoes produces the best estimate of 
the lightning locations for subsequent analysis. These lightning locations are next clus- 
tered into groups that define single thunderstorm cells and storm complexes. This is ac- 
complished by creating grids of lightning activity in 5-10 min intervals and finding iso- 
lated lightning density maxima within the grids. These maxima, typically 2 flashes per 
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100 km 3 , are then used to determine the total number of storms and as seed (initial 
guess) points for a K- Means nearest neighbor clustering algorithm that optimally assigns 
the flashes to the proper storm. Subsequent groupings of the individual storms permitted 
the identification of entire storm complexes. Uncorrected location errors will degrade 
the process of objectively identifying and assigning the individual lightning strikes to 
their parent thunderstorms These misclassifications, in turn, could alter the true nature 
of the lightning life-cycle time series. 

Lightning data offers synergistic information to the radar tracking and cell iden- 
tification algorithms in operational use. The NEXRAD radar storm identification tech- 
niques rely on a base scan reflectivity threshold to identify storm cells. It has been 
demonstrated that lightning strikes can be clustered into discrete cells, thus offering an 
opportunity for augmenting the decision criteria that are based on the radar algorithms 
alone. 
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CHAPTER VII. 


RECOMMENDATIONS 

During the 1990s, the NEXRAD radars will provide storm coverage over the 
conterminous United States (Figure 70). Figure 71 shows the radar coverage at a range 
of 125 km with the track of the 15 November 1989 Huntsville, AL tornado superim- 
posed. At 125 km, for example, the 1° NEXRAD radar beam is about 1 km in diameter 
and the base-scan height (minimum elevation) is 1 km above the surface of the earth. 
Note that the distance from the Nashville, TN and Birmingham, AL NEXRAD radars to 
the tornado is 190 km. At this distance the beam width is 3.4 km. Such a broad sam- 
pling volume will degrade the performance of the storm identification and tracking al- 
gorithms, and the Doppler wind velocity estimates will be compromised. Yet, the 
ground-based lightning network will provide overlapping coverage with the radar for 
improved storm identification. 

Finally, a lightning mapping sensor using CCD-focal plane technology will allow 
total lightning activity to be monitored continuously from space (Christian et al., 1989). 
This instrument is planned for the next generation of operational weather satellites 
(called GOES-Next) and will provide coverage from Canada to Brazil (Figure 72). These 
data will have a spatial resolution of 10 km, sufficient for the identification of in- 
dividual thunderstorms. The satellite lightning data will not require the sensitive site er- 
ror corrections required by the cloud-to-ground lightning networks. The clustering and 
extrapolation forecasting methodologies are equally applicable to the data collected by 
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Figure 70. NEXRAD network sites. 
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NEXRRD COVERAGE AT 125 KM 



Figure 71. NEXRAD sites in the Southeastern United States with 125 km range 
circles. The track of the 15 November 1989 tornado at Huntsville, Alabama 
is indicated by the cross and solid line. 
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Figure 72. Proposed 10.5° field of view centerd at 2°N latitude and 75° longitude 
for the lightning mapper sensor on GOES-Next. 
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the lightning mapper. Sensor fusion with lightning data should provide more capability 
for diagnosing the present and future weather situation than any one of these sensors can 
possibly offer by itself. 
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APPENDIX A: THE FFIX ALGORITHM 


Mathematical Formulation of the Best Poi nt Estimate (BPE1 

FFIX is based on the result that if there are no bearing errors, then the target 
vector T lies in each of the individual bearing planes and, therefore, is normal to the 
normal vectors, N., to each of the bearing planes, i.e., TN.-O (Figure 73). The ob- 
served bearing planes are defined by their respective bearing and station vectors as 
measured from the center of the earth. The bearing vector, 0,, gives the direction from 
which the signal from the lightning strike is sensed at station S.. Due to the random and 
systematic errors, the observed bearing planes do not generally contain a common target 
vector, T. The true bearing to the target is represented by p r Thus, we choose T such 
that the dot products, T N { , are a minimum. FFIX applies the method of weighted least 
squares that minimizes the objective function 

f(T) « w.(T'N/, (A.l) 

subject to the constraint |T| ■ 1. The choice of the weights, w., is based on the prin- 
ciple of maximum likelihood first examined by Stansfield (1947). This is equivalent to 
minimizing 

s k “ E < £ i /£ f (A.2) 
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where is the sum of squared bearing errors for a set of k bearings, fj is the difference 

* . 

between the observed bearing and the bearing from station S ; to the BPE and a is the 
range weighted bearing standard deviation for the ith station. That is, 

a - a.Dj, (A.3) 

where cr is the standard deviation of the bearing errors and D. is the distance from DF 
station S ; to the BPE. For the first iteration when there is no BPE the weight for DF 
station Sj is simply 

w. = (l/<7j) 2 . (A. 4) 

Each successive iteration uses the previous BPE for its initial location estimate, 
which in turn determines each of the D.s. Thus, DF stations having larger bearing er- 
rors and at a further distance from the BPE will be given lesser weight. Define 

T as a 3 x 1 column vector, 

N as a n x 3 matrix of normals 
W as a n x n diagonal matrix of weights. 

If C is defined as C « N T WN, then the objective function (Eqn. A.l) can be expressed 
as 

f(T) =* T t CT . (A. 5) 

For real data C will be positive definite, i.e., f(T) > 0 if T ^ 0. 
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Simplifying the Problem 


First, T t CT - k is the equation of an ellipsoid for any k > 0. The problem is 
simplified if we rotate the axes of the ellipsoid to get the equation in standard form (see 
also Gething, 1978). Note that 

T - Py (A.6) 


is a rotation if 


P T P - I, 


(A.7) 


where P is the matrix of rotation and I is the identity matrix. The required rotation 
gives 


T t CT - y T P T C Py = y T A y, (A.8) 


where 


Thus, 


A - diagonal (A p A r A s ) . 


f(T) - A x y* + A,^ + A 3 y^ , 


0 < Aj <_ A 2 <_ Aj 


The minimum value of f(T) is A x for y * (1,0,0) or 


(A.9) 


(A.10) 


T - Py 


£ 

r 31 


(A.ll) 


The problem is reduced to determining A x and P. The multipliers A. must satisfy 
the equation 
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|C - AI| -0 . 


(A. 12) 


Eqn. (A. 12) is the characteristic polynomial of C and the minimum root A x can be found 
directly without an iteration procedure. The solution vector T can now be computed 

from the matrix equations C « AI and | T | - 1. 

From the geometry in Figure 73 and the properties of the dot and cross products 

we get 


(A. 13) 
(A.14) 

For the spherical triangle defined by the points TSjN. we get the relationship 


T'Nj - cosoj 
| T S. | - sin*. . 


Therefore, 


cosa. « sinpj sine. . 


(A. 15) 


(T'Nj) 2 - cos J Oj 

- | T x S. | 2 sin 2 <r. . 


(A. 16) 


From this result it follows that by defining 

Wj - (1/<7* 2 )|TxS.J 2 (A.17) 


the function 

E(sin 2 t-Jo* 2 ) =* Ht-Jo') 1 (A. 18) 

is minimized. The error introduced by the approximation sincj - e. is < 1% for ^ * 14°. 
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The Confidence Region 


The confidence ellipse calculation is based on the perturbations (variance) of the 
bearings to the BPE, i.e., bearing errors and target position errors are treated as dif- 
ferentials (see Stansfield, 1947). Linearization of the functional relationship between the 
BPE and observed bearings gives an approximate expression for the inverse covariance 
matrix of T. Therefore, the bearings to the BPE are used, not the observed bearings. 
Since both T and Sj lie in the bearing plane, the normal vector can be calculated from 

N. - (T x S;)/|T x Sj | . (A.19) 

As before, C * N T WN, with the exception that now C has the value of zero for 
an eigenvalue. Since T'N. = 0, or in matrix terms, NT = 0 implying that CT * 0, the 
equation 

x T Cx - k 2 (A.20) 

can be rotated by x = Py to 

y T A y - A 2 y 2 2 + A 3 y 3 2 - k 2 (A.21) 

where P T CP = A - diagonal (0,A 2 ,A 3 ). 

When the site errors have been removed, the random errors are normally dis- 
tributed and the errors with respect to the lines of position (i.e., the BPE) are independ- 
ent. Eqn. (A.21) represents an ellipse in the y 2 y 3 plane which is the plane tangent to the 
earth at T. If x is the deviation from T, then x is normally distributed in a plane tan- 
gent at T, and x Cx has the chi-square distribution with 2 degrees of freedom. 

Given the probability (x 2 <_k 2 ) * P for 

k 2 = -2 In (1-P), (A. 22) 
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the axes of the confidence ellipse for probability level P are given by 


a - (Rk/^A 2 ) ; b - (Rk/vA,) (A.23) 

where a and b are the semi-major and semi-minor axes of the ellipse in units of km and 
R - 6378 km is the radius of the earth. The major axis lies along the eigenvector e 2 
corresponding to A 2 , and the orthographic projection of e 2 onto the plane tangent to the 
earth at T determines the orientation of the confidence ellipse. If T ■ (t p t 2 , t 3 ) and e 2 - 
(a r a 2 , a 3 ) and 8 is the bearing of the major axis, then 

sin8 - (tja 2 - tja^/v^tj 2 + t 2 2 . (A.24) 

The ellipse represents a region of minimum area associated with a given proba- 
bility that the lightning strike occurred on or within the ellipse perimeter. The values of 
k 2 for the 90% (P-0.9) and 50% (P-0.5) ellipses, respectively, are 4.61 and 1.39. (See 
Gething, 1978 for more discussion on the confidence ellipse.) 
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APPENDIX B: LIGHTNING ANALYSIS SOFTWARE 


The lightning analysis software consists of programs to decode the hexadecimal 
data archived on 1600 bpi magnetic tape, convert the data into geophysical quantities, 
perform data quality control, remove site errors, compute the optimal flash position, and 
compute error statistics. 
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c 


SUBROUTINE MAI NO 


C 

C LLP LIGHTNING DATA DECODER, USES THE DDAH TAPE FORMAT C 
C CONVERTS RAW DATA TO PROPER UNITS , REMOVES SITE C 
C ERRORS FOR THE MSFC NETWORK, CALLS FFIX FOR OPTIMAL C 
C LOCATION ESTIMATES WHEN 2 OR MORE DPS DETECT THE C 
C FLASH, CALLS NBTFIX FOR TWO BEARING CUT SOLN C 
C C 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DIMENSION IDATA(128) 

CHARACTER G DATE *8 
CHARACTER* 1 CR , LF , AT 

CHARACTER* 1 ICHAR1 (512) , ICHAR2 (46) , ITMP (46) 

CHARACTER* 2 NBUF(256) 

EQUIVALENCE (IBUF, ICHAR1) , (ITMP(l) ,ICHAR2) , (NBUF, IBUF) 
EQUIVALENCE (IBUF,IBUF1) 

COMMON/ BUFS/ IBUF I (128) , IBUF2 ( 128 ) 

COMMON/ COM1/ 1 BFLAG,NREC, INCHX, IN (5) , IN2 (5) 

COMMON/ SO LN/SNKA, SSGNL, TLAT , TLON , RADIUS , AREA, TSOLN , 

* IYDSN , I YR , MON , NDAY , NRS , NBR , IRES , NEVT , SMA , SMI , ORIEN 
COMMON/ DFSTUF/DFE1 (128 ) ,NS1(128) ,DFR10(128) ,NS10(128) , 

*DFR15 ( 128) , NS15 (128) ,SQ1(128) ,SQ10(128) ,SQ15(128) 
COMMON/SUMRY/KTWOER,NTHRS(4) , NBADF , NBDF2 , NDUP ( 4 ) ,NOVR(4) , 
*NDFHT{4) , KDVG 

COMMON/ KOUNTR/KNT1 , KNT2 , KNT3 , KNT4 , KNTHH , MTKNT , IDFTST, ISECOD 
COMMON/ ECNST/TMX , RSML1 , RSML2 , RBIG, SBIG, BINSIZ 
DATA CR, LF, AT/ZOD, Z25 , Z7C/ 

DATA KLl,KL2,ISTAT,MMNUM/26, 46, 0,8000/ 

CALL DATEG (GDATE ) 

CALL TIME ( ITIME) 

IHTIM- INT { ITIME/ 3 60 000 ) 

AHTIM- FLO AT ( ITIME ) / 3 6 0 0 0 0 . 

AMTIM- (AHTIK-IHTIM) *60 . 

IMTIM-INT(AMTIM) 

AS TIM— 60 . * (AMTIM-IKTIM) 

STIM-FLOAT((IHTIM*100 +IMTIM)*100) + AS TIM 
WRITE (6,3) GDATE, STIM 

3 FORMAT (1H1/20X, 'LLP QUALITY CONTROL ANALYSIS '//20X, '(' ,A8, 2X, 

*F10 . 2 , ' ) '/) 

C 

C THESE VARIABLES ARE USED AS CONSTRAINTS IN ERSTAT FOR SITE ERRORS 
IF NOT WANTING SITE ERRORS, SET IDFTST-0, AND ISECOD-2 


IDFTST-0 
TKX-20 . 

RSML1-30. 

RSML2-50. 

RBIG— 200 • 

SBIG— 5 . 

BINSIZ— 6 . 

C 

C SITE ERROR CODES: O-NONE, 1-SITERC, 2-SITER2 , 3-SITERR 
ISECOD— 2 
IOUNIT— 10 
C 

DO 56 L-1,128 
DFE1 (L) -0 . 

NS1 (L) —0 
DFRIO(L) -0. 

NS10 (L) —0 
DFR15 (L) —0 . 

NS15 (L) —0 . 

SQ1 (L) —0 . 

SQ10 (L) —0 . 

SQ15 (L) —0. 

56 CONTINUE 

DO 66 L-1,5 
IN2 (L) —0 
66 IN (L) — 0 

C 

C.. INITIALIZE COUNTERS 
C 

KNT1-0 

XNT2-0 

KNT3-0 

KNT4-0 

NBADF— 0 

NEVT-0 

KDVG— 0 

KTWOER-O 

MTKNT-0 
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DO 68 L— 1 , 4 
C NQC(L)-0 

NDFHT<L)-0 
NTHRS(L) -0 
68 CONTINUE 

IBFLAG-0 
LEN- 5 12 
INCHK-0 

5 IF(ISTAT.EQ.l) CALL ERSUH 

C DO 5 IJK-1,MMNUM 

C WRITE (6,72) IJK 

C72 FORMAT (IX/ 2 OX, 'RECORD NO.* ',16/) 

C IF(ISTAT.EQ.l) CALL ERSUM 

CALL GETREC (LEN, ISTAT) 

NSIZ-0 

NEXT-1 

DO 10 1-1 , LEN 
C I DATA ( I ) -IBUF (I) 

C IPT— 2*1 

IPT— I 

NSIZ-I-NEXT+1 

C WRITE (6, 11) I , NSIZ , NEXT, ICHAR1 (I) ,ICHAR1(I) 

11 FORMAT (IX, ' I ,NSIZ , NEXT, I CHARI. • ' , 314 , A2 , 2X, Z2) 

IF ( ICHAR1 (I) . EQ.LF. AND. NSIZ. EQ.KL1. AND. IBFLAG. EQ.O) THEN 
CALL DFDATA (ICHAR1, IPT, IMES) 

NEXT— 1+1 

IF(IMES.EQ.l) GOTO 10 

C ELSE IF(ICHAR1 (NEXT) . EQ. AT.AND.NSIZ.EQ.KL2)THEN 

C CALL PADATA(ICHAR1,IPT1 


C 

49 


C 

400 


C 

c 


10 

c 

c 

c 

c 

c 


900 


50 

C 

500 


NEXT— 1+1 

IF ( I CHARI ( I ) . EQ . LF . AND . NS I Z . NE . KL1 . AND . I BFLAG . EQ . 0 ) THEN 
WRITE (6 , 49) I ,NSIZ , NEXT, ICHAR1 ( I) , I CHARI (I) 

FORMAT (IX, ' . . . . I , NSIZ , NEXT , ICHAR1 . ,',3I3,A2,2X,Z2) 

NEXT— 1+1 
GO TO 10 
ENDIF 

SORTING AND SAVING DATA FROM NEXT BUFFER TO COMPLETE 
A PARTIAL STRING 

IF (I BFLAG . NE . 0 ) THEN 
I BFLAG- I BFLAG+ 1 
ISAVPT— 2*1 BFLAG 
ISAVPT-IBFLAG 

WRITE^ 6^400^ C I , IPT^NSIz! IBFIAG, ICHAR2 (IBFIAG) , ICHAR2 (IBFIAG) 
FORMAT ( IX, '** I , IPT, NSIZ, IB, ICHAR2 () ' ,4I4,A2,2X,Z2) 

ENDIF 

IF (ICHAR2 (IBFLAG) . EQ . LF . AND . I BFLAG . EQ . KL1 ) THEN 
CALL DFDATA { I CHAR2, ISAVPT, IMES) 

NEXT— 1+1 

IF (IMES . EQ. 1) GOTO 10 
IBFLAG-0 

ELSE I F ( I BFLAG . EQ . KL2 ) THEN 

ELSE IF (ICHAR2 (IBFLAG) . EQ. AT. AND. IBFLAG. EQ. KL2) THEN 
CALL P ADATA( I CHAR2, ISAVPT) 

NEXT— 1+1 
IBFLAG-0 
ENDIF 
CONTINUE 

SAVE THE LAST DF OR PA PARTIAL DATA BLOCK FROM BUFFER I 
TO COMPLETE THE CHARACTER STRING WHICH CONTINUES IN THE 
NEXT BUFFER 

IF(ICHAR1 (LEN) . NE . LF ) THEN 
K— 1 

DO 50 J-NEXT, LEN 
I F ( K . GT . KL2 ) THEN 
WRITE ( 6 900) 

FORMAT (//IX, 'NO LINE FEED AT END OF RECORD') 

GOTO 5 
ENDIF 

ITMP(K) — ICHAR1 ( J) 

WRITE (6, 500) NEXT, LEN, I BFLAG, K,ITMP(1) ,ITMP(K-1) 

FORMAT ( IX , ' ! 1 NEXT, LEN, IB, K, ITMPK, ITMPNXT ,4I4,2A2) 

I BFLAG- LEN-NEXT+ 1 
I SAVPT-2* IBFLAG 
ISAVPT-IBFLAG 
ENDIF 
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C5 CONTINUE 

GOTO 5 

C CALL ERSUN 

C STOP 

END 

SUBROUTINE GETREC (LEN, ISTAT) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


C c 
C RETRIEVES A 512 BYTE (CHARACTER) RECORD FROM DISK C 
C AND RETURNS TO MAINO C 

c c 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

cc 

IMPLICIT DOUBLE PRECISION (A-H, O-Z) 

EQUIVALENCE (IBUF,IBUF1) 

COMMON/BUFS/IBUFl ( 128) , IBUF2 (128) 

CHARACTER* 20 DFFIL 

INTEGER BFLN 

BFLN-128 

NREC-0 

IBFLAG-0 

READ(9, 2, END-400) (IBUFl(IO) ,10-1,128) 

2 FORMAT (128A4) 

NEXT-1 
LAST- LEN 

C WRITE (6 , 300) (IBUF1 (10) ,10-1, 128) 

300 FORMAT (IX, 3 2 A4) 

GOTO 500 
400 ISTAT- 1 

500 RETURN 
END 

SUBROUTINE DFDATA (LINE, IPT) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C c 

C UNPACKS AND DECODES RAW DF DATA ( DDAH FORMAT) FROM TAPE C 
C C 

C 

IMPLICIT DOUBLE PRECISION (A-H, O-Z) 

C LINE CONTAINS RAW DF DATA STRING W/ CRLF 

C TMP CONTAINS SUBSTRINGS FOR SWAPPING MSB/ LSB 

CHARACTER*! LINE (512) 

CHARACTER TMP* 4 
CHARACTER* 1 OUTBUF(24) 

CHARACTER* 2 4 IOBUF 

CHARACTER* 1 I A POL 

COMMON/ BUFS/ 1 BUF1 (128) , IBUF2 (128) 

COMMON/ COM1/ I BFLAG , NREC , INCHX , IN ( 5 ) , IN2 ( 5 ) 

COMMON/ DFSTUF/DFE1 (128) , NS1 (128) , DFR10 (128) ,NS10(128) , 

*DFR15( 128) , NS15 ( 128 ) ,SQ1(128) ,SQ10(128) ,SQ15(128) 

COMMON/ SOUf/SNKA , SSGNL , TLAT , TLON , RADIUS , AREA , TSOLN , 

* I YDSN , I YR , MON , NDAY , NRS , NBR , IRES , NEVT, SMA, SMI , ORIEN 
COMMON/ SUMRY/KTWOER , NTHRS ( 4 ) , NBADF , NBDF2 , NDUP ( 4 ) ,N0VR(4) , 
*NDFHT(4 ) , KDVG 

COMMON/ KOUNTR/ KNT 1 , KNT2 , KNT3 , KNT4 , KNTHH, MTKNT , IDFTST, ISECOD 

COMMON/ ECNST/TMX , RSML1 , RSML2 , RBIG , SBIG , BINSIZ 

INTEGER IODF{24) , DF1TST 

INTEGER HX80 

INTEGER*4 JD , JDAY , I YEAR 

EQUIVALENCE (IBUF,IBUF1) 

EQUIVALENCE (OUTBUF, IOBUF) , (NCNT,NBR) , (IERR, IRES) 

DIMENSION IRS (4) , IHH (4) , IKM(4) ,ISS(4) ,IMS(4) , ATOP ( 4 ) , IYYDDD ( 4 ) 
DIMENSION POL( 4 ) , CBRG2 (4) ,A2M(4) ,AMP(4) ,TM(4) ,CBRG1(4) ,AZM2(4) 
DIMENSION ITRS (4 ) , TAMP (4) ,B(4) ,IA(2) , AB(2) 

DIMENSION DLATI ( 4 ) , D LONG I ( 4 ) 

DATA DLATI/34. 649167, 35. 399167, 35. 83750, 34. 716667/ 

DATA DLONGI/86. 669167, 86. 076944, 87. 443889, 87. 881667/ 

DATA 1ST, MNODF/26 , 4/ 

DATA HX80/Z80/ 

DTR-O. 01745329 
EMISS1— 999.99 
EBAD1— 888.88 
ERAD1— 6371 . 

ISTART— IPT-IST+1 
IEND-IPT 

C WRITE ( 6 , 200) ISTART, IPT, (LINE (IS) , IS— ISTART, X END) 

200 F0RMAT(1X, 'ISTART, IPT', 214, 2X,26A4) 

II-l 

J— ISTART 
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c 

c .. 
c 


WRITE DF STRING INTO OUTBUF FOR INTERNAL I/O, SWAPPING HSB, LSB 


OUTBUF(l) *LINE ( J) 

OUTBUF ( 2 ) -LINE ( J+ 1 ) 

OUTBUF { 3 ) -LINE (J+4 ) 

OUTBUF ( 4 ) -LINE ( J+5) 

OUTBUF ( 5 ) -LINE ( J+2 ) 

OUTBUF ( 6 ) -LINE ( J+3 ) 

OUTBUF ( 7 ) -LINE ( J+8 ) 

OUTBUF ( 8 ) -LINE ( J+9 ) 

OUTBUF ( 9 ) -LINE ( J+ 6 ) 

OUTBUF ( 10 ) -LINE ( J+7 ) 

OUTBUF (11) -LINE ( J+12) 

OUTBUF ( 12 ) -LINE (J+13 ) 

OUTBUF ( 13 ) -LINE ( J+10) 

OUTBUF (14) -LINE (J+ll) 

OUTBUF ( 15) -LINE ( J+14 ) 

OUTBUF ( 16) -LINE ( J+15 ) 

OUTBUF ( 17 ) -LINE ( J+ 18 ) 

OUTBUF ( 18 ) -LINE ( J+ 19 ) 

OUTBUF ( 19 ) -LINE ( J+ 16 ) 

OUTBUF ( 20 ) -LINE (J+17 ) 

OUTBUF (21) —LINE ( J+22) 

OUTBUF ( 2 2 ) -LINE (J+2 3 ) 

OUTBUF (23) -LINE (J+2 0 ) 

OUTBUF (24) -LINE (J+2 1) 

C 

C.. CHECK FOR ILLEGAL CHARACTERS IN DF DATA STRING 
C 

IMES-0 

DO 30 KC-1,24 

I F ( OUTBUF ( KC) * LT. 'A' .OR. OUTBUF (KC) .GT. ' 9 ' ) THEN 
WRITE (6, 455) OUTBUF 

455 FORMAT (5X, ' BAD CHARACTER IN DF STRING - ',24A1) 

C WRITE (6,230) ( IBUF1 (10) , 10- 1, 128) 

C230 FORMAT (IX, 3 2 A4) 

IMES-1 

RETURN 


30 

191 

C 

C 

C. . 
C 


553 


ENDIF 

CONTINUE 

READ (UNIT-IOBUF , FKT- 191) 101,102 , 103, 104 , 105 , 106, 107,108 
FORMAT (Z2 , 24 , Z4 , Z4 , Z2 , Z4 , Z2 , Z2 ) 

WRITE (6,193) 101 , 102 , 103 , 104 , 105 , 106 , 107 , 108 

GET DF IDENTITY, IDF-I01+1 


IDF— 101+1 

IF (IDF. LT. 1 . OR. IDF.GT.MNODF) THEN 
IQCKEY— 1 

WRITE (6, 553) OUTBUF 

FORMAT (5X, ' ILLEGAL DF ID NUMBER IN STRING - ',24A1) 
RETURN 
ENDIF 


C 

c. . COUNT SUM OF EACH DF DETECTIONS 
C 

IF ( IDF. LE . MHODF) NDFHT (IDF) — NDFHT ( IDF) +1 

C. . GET FLASH POLARITY HEX80 ADDED TO BYTE 12 ON TAPE IF POSITIVE 
C 

POL(IDF)— 1.0 

IF (107 .GE.HX80) POL(IDF)-1.0 

C. . GET AMPLITUDE, 1500 IS OVERANGE VALUE (DF SIGNAL SATURATION) 

^ AMP (IDF) —POL (IDF) * (M0D(I07 , HX80) *256. + FLOAT ( I08))/10. 

IF (ABS (AMP (IDF) ) . GE. 1500 . ) THEN 
NOVR ( IDF ) — NOVR ( IDF ) +1 
C 

C WRITE (6, 143) IDF, AMP(IDF) 

C143 FORMAT ( IX , ' DF OVERRANGE' , 14 , F8 . 2) 

RETURN 

ENDIF 

IF (ABS (AMP (IDF) ) . LT . 10 . ) THEN 
NTHRS(IDF) - NTHRS (IDF) +1 
RETURN 
ENDIF 

C.. GET MINUTE OF THE DAY AND CONVERT TO HHMM FORMAT 
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IHH ( IDF) -INT( 102/60) 

IMM ( IDF) “102*60 *IHH( IDF) 


151 





c 

c. . 
c 


c 

c.. 

c 


c 

c. . 
c 


300 


350 


400 

450 

C 

c. . 
c 

c 

c.. 

c 

c 

c. . 
c 

c 

c. . 
c 

c 

193 

c 

C. . 

c 


55 


C 

c. . 
c 

c 


196 


230 

C 


GET MILLISECONDS OP THE MINUTE 

ISS (IDF) —INT( 103/1000) 

IMS ( IDF) *103-1000* ISS ( IDF) 

GET TIME OF DAY HHMMSS FORMAT 

TM(IDF) -FLOAT ( ( (IHH(IDF) *100) +IMM(IDF) ) *100+ISS (IDF) ) 
*+FLOAT(IMS (IDF) )/1000. D00 ' U 

GET THE YEAR 


JDAY-I04+1 

NYEAR-JOAY/365 

LP YEAR- NYEAR/ 4 

J D-LPYEAR+ 3 6 5 *NYEAR 

IF ( JOAY-JD) 350,300, 400 

JDAY-365 

KYEAR— LPYEAR*4 

IF (XYEAR. EQ . MYEAR) JDAY-366 

NY EAR— NY EAR- 1 

GOTO 450 

NYEAR— NYEAR- 1 

LP YEAR- NY EAR/ 4 

JCHLPYEAR+ 3 6 5 * NYEAR 

JDAY— JDAY-JD 

IF(JDAY.EQ.O) GOTO 300 

GOTO 450 

JDAY-JDAY-JD 

JDAY-JDAY-1 

GET YYDDU FORMAT FROM NYEAR AND JDAY 
lYYDDD(IDF) -NYEAR* 1 00 0+J DAY 
GET THE EQUIVALENT MONTH, DAY, AND YEAR 


GET NUMBER RETURN STROKES 


IRS (IDF) -105 

GET UNCORRECTED DF AZIMUTH TO FLASH 
AZM(IDF) — (FI0AT(I06)/65536. ) *360. 

FORMAT ( IX , 'IODF ',818) 

CHECK DF DATA FOR UNACCEPTABLE VALUES 
IQCKEY-0 

IF (IDF. LT. 1 . OR. IDF. GT .MNODF) THEN 
IQCKEY— 1 

WRITE (6 , 55) OUTBUF 

FORMAT (5X,' BAD CHARACTER IN DF STRING - ' 24A1) 

RETURN 

ENDIF 

IF ( IHH(IDF) .LT. O.OR. IHH(IDF) .GT.24) IQCKEY— 1 

IF(IMM(IDF) .LT. O.OR. IMM( IDF) .GT.60) IQCKEY— 1 

IF (ISS (IDF) .LT. O.OR. ISS (IDF) .GT.60) IQCKEY— 1 

IF(IMS (IDF) .LT. O.OR. IMS (IDF) .GT.1000) IQCKEY— 1 

IF ( IYYDDD ( IDF) .LT. 83000. OR. IYYDDD(IDF) .GT. 90000) IQCKEY— 1 

IF ( IRS (IDF) . LT . O.OR. IRS ( IDF) .GT.14) IQCKEY— 1 

IF (AZM(IDF) . LT . 0 . .OR.AZM(IDF) .GT.360.) IQCKEY— 1 

IF DATA UNACCEPTABLE DO NOT CONTINUE WITH STRING PROCESSING 

IF ( IQCKEY . EQ .- 1 ) THEN 
NQC ( IDF) -NQC ( IDF) + 1 

WRITE (6,196) IQCKEY , IDF , IHH ( IDF) , IMM ( IDF) , IS S ( IDF) , IMS ( IDF) . 
► TM(IDF) , IYYDDD(IDF) , IRS (IDF) , AZM(IDF) , POL(IDF) ,AMP(IDF) 
FORMAT ( IX, ' $196 DF DATA ',12,' DF- ',12 , 2X, 312 , 13 , 4X,F10.3, 

* I6,I4,F10.4,F«.1,F8.2/) 

WRITE (6,230) (IBUF1 (10) , IO-l, 128) 

FORMAT (IX, 3 2 A4) 

STOP 

RETURN 

ENDIF 
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c 

c. . 
c 


188 


C 

C 

c 

c 

c. 

c 


c 

C211 


C 

C 

C199 

C 

271 

299 

C 

C. . 

C 

285 

C 

C. . 
C. . 
C. . 

c. . 
c 


c 

c. 

c 

c 


707 


711 

C 

C. . 
C 


GET A CORRECTED BEARING ANGLE FOR THIS DF 
IF {AZM ( IDF) .EQ.O. .AND. IDFTST.EQ. IDF) THEN 

WRITE ( 6 # 188 ) IDF,AZM(IDF),TM(IDF),AMP(IDF) 

FORMAT { IX, ' ZERO DEGR. ',I2,3F12.3) 

CBRG1 (IDF) "0 . 0 
GOTO 285 
ENDIF 

CBRG1 ( I DF) -AZM ( IDF) +S ITERC ( IDF , AZM ( IDF) ) 

CBRG1 (IDF) —AZM (IDF) +SITERR ( IDF, AZM ( IDF) ) 

CBRG1 (IDF) — AZM(IDF)+SITER2 (IDF, A2M(IDF) ) 

CBRG 1 ( I DF ) - AZM ( I DF ) 

CHECK THAT CORRECTED ANGLE NON NEGATIVE, ADD 360 FOR PROPER ONE 

IF(CBRG1 (IDF) . LT.O.) THEN 

WRITE (6,211) CBRG1(IDF),IDF,AZK(IDF) 

FORMAT (IX, 'ANGLE LESS THAN 0. . , F12 . 3 , 14 , F12 . 3/) 

CBRG1 ( IDF) — CBRG1 (IDF) +360 . 

FNDIF 

IF (TM(IDF) .GT. 225300. .AND.TM(IDF) . LT. 225500. ) THEN 
WRITE (6, 199) TM(IDF) , IDF, AZM (IDF) ,CBRG1(IDF) 

FORMAT ( IX, ' DF, BEARINGS ' , F12 . 3 , 2X, 12 , 2F10 . 4 ) 

ENDIF 
GOTO 285 

FORMAT (IX, 'NO CORRECTION YET AVAILABLE FOR DF1') 

CHECK FOR NUMBER OF DFS IN TIME COINCIDENCE 

CALL TIMECO (IDF,TM(IDF) , ICNT1, ICC, IYYDDD) 

FIND NUMBER OF DFS IN COING^ENGE.^FJDFS-^ CA^ NETFIX 
IF # DFS IS 2 OR MORE USE OPTIMIZER FFIX. IP FFIX FAILS 
TO GET A GOOD SOLN, THEN GET BEST CUT FROM NETFIX BASED 
ON MIN. SEMIMAJOR AXIS OF THE ERROR ELLIPSE. 

IF ( IN ( 5) . NE . 1) GOTO 275 
IERR-0 

IF ICNT1 IS GREATER THAN 4, ONE OR MORE DFS SAW MORE THAN ONE 
EVENT WITHIN THE TIME COINCIDENCE WINDOW 

IF ( ICNT1 .GT . 4 ) THEN 

StU^'Df'sM 1 ^ 1 ™^ ONE FLASH WITHIN TIME WINDOW. 
* 14 , F12 .3) 

ICNT1-4 
ENDIF 
NCNT— ICNT1 
ISUM-0 

DO 711 IK-1, NCNT 
DO 711 IK-1, MNODF 
IF ( IN ( IK) . EQ . 1 ) THEN 
IN2 ( IK) —IK 
ISUM—ISUM+1 
AZM2 ( IK) —AZM ( IK) 

CBRG2 (IK) — ATMP ( IK) 

ENDIF 

CONTINUE 

CHECK TO SEE IF NUMBER OF DFS IN SOLN EQUALS COUNTER NCNT 


7277 


C 

c 

c. 

c 


c 

2112 


1I(I1 „„ 

FORMAT (/IX, 'TIME CORR ERROR AT , F12 . 3 , 12 (12, IX) ) 
GOTO 901 

“°TM ( IDF) . GT . 185500 .. AND . NDAY . EQ . 14 ) STOP 

USE THIS FOR COMPUTING DF SITE ERRORS 

IF(IDFTST.EQ.O) GOTO 75 
IF(IN2 (IDFTST) .EQ.O) GOTO 901 

150 IF (IN2 (I) .EQ. IDFTST. AND. NCNT. GE. 3) THEN 
IN2 (I) —0 

TAZ1-CBRG2 (IDFTST) 

NCNT-NCNT-1 
ICNT1— NCNT 
ENDIF 


0 
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717 CONTINUE 

IF (NCNT.EQ.2 . AND.IN2 (1) .EQ.IDFTST) GOTO 901 
IF(NCNT.EQ.2 . AND. IN2 (2) .EQ.IDFTST) GOTO 901 
C 

75 DO 712 KI-l,MNODF-l 

IF (IN2 (KI) .EQ. 0. AND. IN2 (KI+1) .NE.OJTHEN 
IN2 (KI) — IN2 (KI+1) 

A2H2 ( KI ) -AZM2 ( KI+ 1 ) 

CBRG2(KI)-CBRG2 (KI+1) 

IN2(KI+1)-Q 

C IF (TM(IDF) .GT. 185200 . .AND.TM(IDF) .LT. 165500. ) THEN 

C WRITE (6,76) KI , NCNT , IN , IN2 , AZM2 , CBRG2 

76 FORMAT (IX, 'DEBUG KI ',12 (12 , IX) , 8 (F7 . 3 , IX) ) 

C ENDIF 

ENDIF 

712 CONTINUE 

C IF (TM(IDF) .GT. 201815. .AND.TK(IDF) . LT. 201818 . )THEN 

C WRITE (6,76) KI , NCNT , IN , IN2 , AZM2 , CBRG2 

C ENDIF 

C 

C. . MAKE A SECOND PASS THROUGH ABOVE SORTING LOOP IF DFS NOT 


C. . LISTED CONSECUTIVELY 
C 

IF(IN2 (1) .EQ. O.OR.IN2 (2) .EQ. 0) GOTO 75 
C IF (NEVT.GE. 80) STOP 

C 

IF ( ICNT1 . GE . 2 ) CALL FFIX ( IN2 , CBRG2 , TM) 

C IF(ICNT1.GE. 2) THEN 

C811 CALL FFIX ( IN2 , CBRG2 , TM) 

C WRITE (6,848) TSOLN , IERR, TLAT, TLON , SMA, SMI , ORIEN, IN2 , CBRG2 , TM 

848 FORMAT </2X, 'FFIX RETURN ' , F12 . 3 , 12 , 5F12 . 3/1X, 5 (12 , 2X) , 

* 4 (F12 . 4 , IX) , 4 ( F12 . 3 , IX) ) 

C IF (IERR. EQ. -1) GOTO 901 

IF (IERR. NE. 1 . OR. ICNT1 . EQ. 2) THEN 
C WRITE (6,841) TK(4),IERR 

C NBADF-NBADF+1 

C IF(IERR.EQ.O.OR.ICNT1.EQ.2)THEN 

C IF (ICNT1 .GE. 2} THEN 

C 

C.. THIS FORCES BEST 2 DF FIX (BEST - MIN (SMA)) 

C 

IBDCHX-0 

TSTOR3-999. 

SKA-999. 

ICNT1-2 
NCNT— ICNT1 

DO 745 I— 1 , MNODF-1 
DO 754 J— 2 ,MNODF 

IF(I.GE.J) GOTO 754 
IA(1)— IN2 (I) 

IA(2)— IN2 (J) 

IF(IA(1) . EQ. O.OR. IA(2) .EQ.O) GOTO 754 
AB ( 1 ) — CBRG2 ( I ) 

AB(2) — CBRG2 (J) 

CALL FFIX(IA # AB,TM) 

C WRITE (6,872) IER, TSOLN, TLAT, TLON, SMA, IA f AB 

C872 FORMAT (IX, 'LN 872 CHK SOL' , 2X, 12, 4F12 . 3 , 2X, 212 , 2F8 . 2 ) 

IF ( IERR . EQ . 1 . AND . SMA . LE . TSTOR3 ) THEN 
C 

C. . CHECK FOR DIVERGENT ANGLES- SOLN GOES AROUND THE WORLD 
C 

XTL- TLAT*DTR 
XTLO-TLON * DTR 
DO 790 K-1,2 

XSL-DLATI (IA(K) ) *DTR 

XSIO— DLONGI ( IA ( K) ) *DTR 

CALL AZRN (XSL, XSLO,XTL, XTLO, API, RP1) 

ABM1— AB (K) —5 . 

ABP1— AB(K) +5. 

RP— RP1 * ERAD1 
AZP— AP1/DTR 

IF (AZP.GE. ABPl.OR. AZP.LT.ABM1)THEN 
c WRITE (6, 7979) TSOLN, TLAT, TLON, RP, 

c * AZP, IA, AB, SMA 

C797 * FORMAT (IX, 'DVGA' , 4X, 5F12 • 3 , 2X, 

c * 212, 3F8 .2) 

C IF(NEVT.GE.80) STOP 

KDVG-KDVG+1 


790 


GOTO 754 
ENDIF 
CONTINUE 
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754 

745 


775 

C 

C 

c 

870 


C 


641 

C 


C 

C 

C 

c 

C. 

c 

c 


TSTORl-TLAT 

TSTOR2-TLON 

TSTOR3-SKA 

TSTOR4-SMI 

TSTOR5-ORIEH 

TSTOR6-AREA 

TSTOR7 -RADIUS 

IDTl-IA(l) 

IDT2-IA(2) 

TBRl-AB(l) 

TBR2-AB ( 2 ) 

IBDCHK-1 
ENDIF 
CO NT I NUB 
CONTINUE 
TLAT-TSTOR1 
TLON—TSTOR2 
SMA-TSTOR3 
SMI -TS TOR 4 
ORIEN-TSTOR5 
AREA-TSTOR6 
RADIUS-TSTOR7 
IN2(1)~IDT1 
IN2(2)-IDT2 
CBRG2 (1)*TBR1 
CBRG2(2)-TBR 2 
DO 775 KP-3,4 
IN2(KP)-0 
CBRG2 (KP)-O 
CONTINUE 
ENDIF 
ENDIF 

WRITE (6,870) TSOLN, TLAT , TLON , SMA , IN2 , CBRG 2 
FORMAT ( IX , ' IN 870 CHK SOL' , 4F12 . 3 , 2X, 512 , 4F8 .2) 

IF ( IBDCHK . NE . 1 ) THEN 
NBADF-NBADF+1 
GOTO 1801 
GOTO 901 
ENDIF 
ENDIF 

FORMAT ( IX , ' . . • NO FIX FROM FFIX AT',F12.3,' ERROR-', 13) 
IF(IN(1) • EQ. 0) GOTO 901 
XTL-TLAT * DTR 
XTLO— TLON*DTR 
IF ( IDFTST . EQ . 0 ) GOTO 1028 
GOTO 1015 
ENDIF 

IF (ICNT1 . LT. 4 ) GOTO 901 (USE THIS IF ONLY ALL 4 USED) 

SET UP FOR 2 DF FIX 

I F ( ICNT1 . EQ . 3 . AND . IN ( 1 ) . EQ. 1 . OR. ICNT1 . EQ. 2 . AND. IN(1) .EQ.0)THEN 
IF ( ICNT1 . EQ . 2 ) THEN 


C 


715 

510 

C 

C 


C 

1101 

C 

514 

C 

C 

C1801 


WRITE (6,510) ICNT1 , NCNT , IN , IN2 , CBRG 2 , ATMP 
DO 715 Kl-1,4 
IF (ATMP (Kl) . EQ . 0 . ) THEN 
GOTO 715 
ENDIF 

IF ( IN (Kl) . EQ. 1. AND. IN2 (1) . NE . 0 ) THEN 
IN2(2)«K1 
CBRG2 ( 2 ) “ATMP ( Kl ) 

ENDIF 

IF (IN (Kl ) . EQ . 1 .AND. IN2 ( 1) . EQ . 0 ) THEN 
IN2(1)-K1 
CBRG2 ( 1 ) "ATMP ( Kl ) 

ENDIF 

CONTINUE 

ENDIF 

F0RMAT(1X, 'DFDATA. .ICNT1, NCNT, IN, IN2 , CBRG2 , ATMP ' , 12I4/1X, 8F12 . 4) 

IF (ICNT1 . EQ . 2 .AND. IN ( 1) . EQ. 1) GOTO 890 

NCNT-2 

IF(IERR.EQ.O) GOTO 811 
IF(IERR.EQ.-l) GOTO 901 
IF(IERR.EQ.O)THEN 
WRITE (6,1101) NCNT, ICNT1,NBR 

FORMAT (IX, 'HERE WE ARE ',314) 

CALL FFIX ( IN2 , CBRG2 , TM) 

IF(IERR.EQ.-l) GOTO 901 
IF(IERR.NE.l) THEN 


WRITE (6,841) TM(4) , IERR 

CALL NETFIX ( IN2 ( 1 ) , CBRG2 ( 1 ) , IN2 


(2) , CBRG2 (2) , TLAT, TLON) 


TLON— l.*TLON 
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2099 


WRITE(6,2099) TSOLN,IN2 (1) # CBRG2 (1) , IH2 (2) ,CBRG2(2), 
♦TLAT, TLON 

FORMAT (IX, '2 DF FIX' , F12 . 3 ,2X, 12 ,F12 . 3 , 2X, I2,3F12. 3) 
ENDIF 

C IF(IN(1) .EQ.O) GOTO 901 

1801 XTL-TLAT*DTR 

XTLO-TLON*DTR 
IF(IDFTST.GT.O) GOTO 1015 
SKDIF 

IF ( IDFTST . EQ . 0 ) GOTO 1028 
C 

C. . CONVERT SOLN TO RADIANS FOR AZRN TO GET DF1 SITE ERRORS 
C IF DF1 NOT IN SOLN THEN SKIP THIS 

C 

XTIf-TIAT * DTR 
XTLO— TLON * DTR* ( - 1 . 0) 

1015 XSL-DLATI (IDFTST) *DTR 

XSLO-DLONGI ( IDFTST) *DTR 

CALL AZRN(XSL, XSLO,XTL,XTLO, AZ1,R1) 

C 

C.. CONVERT FROM RADIANS BACK TO KM, EARTH RADIUS-6370 KM 
C 

C WRITE (6, 912) AZ1,R1 

912 FORMAT ( IX , 'AZRN SENDS DFDATA AZ1,R1 ',2F10.6) 

A 2 DEG- A Z 1/DTR 
RKM-R1*ERAD1 

C. . GET DF 1 BEARING ERROR 

C SKIP THIS PART IF NOT NEEDING BEARING ERROR 
C 

C GOTO 623 

IF (ATMP( IDFTST) .EQ.O. 00) GOTO 901 
C WRITE (6 ,3111) TSOLN , IN2 , RKM , AZDEG , DIFDEG , CBRG2 

3111 F0RMAT(/1X, 'DIFDEG BEFORE TAZ CHECK ' ,F12.3, 2X, 512, 7F10. 4) 

C TAZ 1-ATMP ( 1 ) 

TAZ 2— AZDEG 

C WRITE(6, 3112) TAZ 1 ,TAZ2 , DIFDEG, XSL, XSLO, XTL, XTLO 

IF(ABS (TA22-TAZ1) . GT. 300. ) THEN 
IF (TAZ 2 * GT . TAZ 1 ) TAZ1-TAZ1+360 . 

IF (TAZ2 . LT. TAZ1) TAZ2-TAZ2+360. 

ENDIF 

DIFDEG-TAZ2-TAZ1 

C WRITE (6 ,3112) TAZ1 , TAZ2 , DIFDEG , XSL, XSLO , XTL, XTLO 

3112 FORMAT (IX, 'AFTER TAZ CHECK ',7F10.4) 

C IF(NEVT.GE.SO) STOP 

C 

C. . BEFORE ACCEPTING THIS SOLN, CHECK FOR LARGE ERROR RADIUS AND 
C LAT/LON SOLN WHICH ARE REASONABLE 
C 
C 

623 IF(SMA.GE. 100 . ) GOTO 625 

IF (TLAT.GE. 42. .OR.TLAT.LT. 30. ) GOTO 625 
IF(TLON.GB. 94. .OR.TLON.LT. 78.) GOTO 625 
625 GOTO 635 

C625 WRITE ( 6 , 680) TSOUf , TLAT, TLON , SMA 

680 FORMAT (IX, 'DISTANT OR BAD SOLN. .REJECT IT AT ',F12.3,2X, 

*2F12 . 6 , ' SEMI -MAJOR AXIS-', F8. 2/) 

C 

C. . CHECK FOR A BAD SOLN USING 2 BEARINGS AND GIVING SAMS 
C SOLN AS PREVIOUS SOI29. THIS IS DUS TO NARROW CUT-ANGLE 
C NBDF2 IS NO. OF WRONG SOLNS FROM FFIX, CUT SOLN BETT ER 
C 

635 IF(ABS (DIFDEG) .GT.TMX) THEN 
NBDF2-NBDF2+1 
GOTO 901 
ENDIF 

I F (TLAT. EQ. PLAT. AND. TLON. EQ.PLON) THEN 
KTWOER-KTWOER+1 

C WRITE (6 , 1777) TLAT , PLAT , TLON , PLON 

C1777 FORMAT ( IX , ' TLAT , PLAT , TLON , PLON- ' , 4 F12 . 6 ) 

GOTO 901 
ENDIF 
C 

C.. SORT DF ANGLE ERROR AS A FCN OF ANGLE IN ERSTAT 
C 

C CALL ERSTAT (AZDEG, RKM, DIFDEG) 

CALL ERSTAT ( TAZ 1 , RKM , DI FDEG , CBRG2 ) 

C 

C.. ESTIMATE THE PEAK CURRENT FOR THIS FLASH WITH FCN FKA 
C 

1028 SKA- FKA (TAMP) 

SNKA-SKA 

C 

C.. FIND THE LARGEST NUMBER OF STROKES IN A FLASH 
C 

156 


ORIGINAL PAGE IS 
OF POOR QUALITY 



oooonoooooooo 


NRS-0 

DO 1411 I-1,NCNT 

1411 I F ( ITRS ( IN2 ( I ) ) . GT . MRS ) NRS-ITRS (IN2 (I) ) 

C 

C NCNT-0 

C 

C.. GET THE EQUIVALENT MONTH, DAY, AND YEAR FOR A FLASH 
C 

CALL YDDMY (IYDSN, NDAY,MON, IYR) 

C 

C.. SET NO RADIUS FROM SOLN TO EMISS1— 999 . 99 IN SOI21 SET 
C 

IF (RADIUS . EQ . 0 . ) THEN 
AREA-EMISS1 
SMA-EMISS1 
SHI-EMISS1 
RADIUS-EMISS1 
ORIEN-EMISS1 
ENDIF 

IF ( RADIUS . GE . 200 . ) THEN 
AREA-EBAD1 
SMA-EBAD1 
SMI-EBAD1 
RADIUS-EBAD1 
0RIEN-EBAD1 
ENDIF 
C 

C.. COUNT NUMBER OF 2 DF AND 3 DF SOLNS 
C 

IF(ICNT1.EQ. 2) KNT2-KNT2+1 
IF(ICNT1. EQ. 3) KNT3-KNT3+1 
IF(ICNT1 . EQ. 4) KNT4-KNT4+1 
IF(ICNTl-EQ.l) KNT1-KNT1+1 
C 

C.. CHECK FOR TIME >24 HOURS, DUE TO DAY CHANGE, NEED TO SUBTRACT 24 
C 

IF (TSOLN.GE. 240000.) TSOLN-TSOLN-24QOOO . 

C.. CORRECT FOR LEAP YEAR, NDAY-NDAY+1 (LLP DAY OF CENTURY COUNT IS 
C.. INCORRECTLY COMPUTING THIS. 

C 

IF (IYR. EQ. 80 ) NDAY-NDAY+1 
C 
C 

C WRITE (6 , 915) TSOLN, IYR,MON,NDAY, AZDEG,RKM, ATKP(l) , DIFDEG, 

C *TLAT , TLON , NRS , SSGNL , SNKA , RADIUS , SMA , SMI , AREA 

C915 FORMAT (IX, ' SOLN AT' , F12 . 3 , 5X, 312 , 2X, 4F12 . 6 , '- DF 1 ERROR'/ 

C *2F12.6,I4,6F12.2) 

C WRITE (6, 918) IN2,CBRG2 

C918 FORMAT ( IX, ' IN2 , CBRG2- ' , 512 , 4 (F6 . 2 , 2X) ) 

C 

C.. WRITE SOLN TO DISK FILE ON EADS (UNIT 10) 

C 

NEVT-NEVT+1 

WRITE (10,3015) NEVT , TSOLN , IYR , MON , NDAY , TLAT , TLON , SNKA , SSGNL , NRS 
3015 FORMAT (I6,2X,F12.3,1X,3I2,2F12.6, 2F10 . 2, 2X, 12} 

WRITE (10,4015) IN2 , CBRG2 , SMA , SMI , ORIEN, RADIUS , AREA 
4015 FORMAT (5I2,4(F6.2,1X) ,5(F10.2,1X) ) 


IF (NDAY. EQ. 14 .AND. TSOLN. GT. 201815 . . AND. TSOLN. LT. 201818 .) THEN 
WRITE (6,230) (IBUF1 (10) , IO-l , 128) 

WRITE (6,915) TSOLN , IYR , MON , NDAY , AZ DEG , RKM , ATMP ( 1 ) , DIFDEG , 

*TLAT , TLON , NRS , SSGNL, SNKA , RADIUS , SMA , SMI , AREA 
915 FORMAT (IX, ' SOLN AT' , F12 . 3 , 5X, 312 , 2X, 4F12 . 6 , DF 1 ERROR'/ 
*2F12.6,I4,6F12.2) 

WRITE (6,3015) NEVT , TSOLN , IYR , MON , NDAY , TLAT , TLON , SNKA , SSGNL, NRS 
WRITE (6 , 4015) IN2 , CBRG2 , SMA, SMI , ORIEN, RADIUS , AREA 
IF (TSOLN.GE. 201818. ) STOP 
ENDIF 

. . WRITE HOURLY SOLN SUMMARY 


IF(NEVT.EQ. 1) THEN 
LSTHH— IHH ( IDF) 

WRITE (6, 2017) 

2017 FORMAT (20X, 'EVT' , 4X, 'TIME' , 5X, ' YYMMDD' , 6X, 'LAT',10X, 'LON', 

* 8X, 'KAMPS' ,5X, 'VOLTS' ,2X, 'RS',4X, 'RADIUS ',4X, 'SMAJ',4X, 

* 'SMIN',3X, 'AREA') 

WRITE (6,3017) NEVT , TSOLN , I YR , MON , NDAY , TLAT, TLON , SNKA , SSGNL, 

* NRS, RADIUS, SMA, SKI, AREA 

3017 FORMAT (/IX, '1ST SOLN (TAPE) , 14 , 2X, F12 . 3 , IX, 312 , 2F12 . 6, 2X, 

* 2F10 .2,I2,4F10. 2/) 

ENDIF 

IF ( IHH ( IDF) .EQ. LSTHH) KNTHH-KNTHH+l 
IF (IHH ( IDF) .NE. LSTHH) THEN 
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LSTH1-IHH (IDF) 

IF ( LSTH1 . LT . 0) LSTHH-0 

WRITE (6, 3020) LSTHH,KNTHH, IYR,MON, NDAY,NEVT 
3020 FORMAT (/IX, 'HOURLY SUMMARY BEGINNING AT ',12, ' UT IS ',I6,4X, 

* 'DATE-', 312,' TOTAL EVENT COUNT SO FAR-', 16 ' 

LSTHH— IHH (IDF) 

KNTHH-0 

ENDIF 

C 

c 

C890 WRITE (6, 899) 

890 CONTINUE 

899 FORMAT ( IX , ' SOLN REQUIRES DF1 FOR 2 DF FIX') 

C 

C. . CLEAR POINTERS 
C 

901 NCNT-0 
ICNT1-0 
LSTDF-IDF 
PREVB— CBRG 1 ( I DF ) 

PREVA- AMP (IDF) 

PREVR-IRS { IDF) 

C IF(LSTDF.NE.l) AZM(l)-0. 

AZDEG— 0. 

RKM-0. 

DIFDEG-0. 

PLAT-TLAT 

PLON-TLON 

TLAT-0. 

TLON-O . 

RADIUS-0 . 

AREA-0. 

SMA-0. 

SMI-0. 

ORIEN-O. 

DO 340 INF-1,5 
IN2 ( INF) -0 
340 IN { INF) — 0 

XTL--1 . 

XTLO— 1 . 

ICC-0 

DO 345 INF-1,4 
B(INF) -0. 

AMP(INF) -0. 

IRS ( INF) -0 
CBRG2 (INF) -0 . 

345 CBRG1 ( INF) —0 • 

CBRG1 ( IDF) -PREVB 
AMP (IDF) -PREVA 
IRS ( IDF) -PREVR 
275 DO 2275 INF-1,4 

ATMP ( INF) — CBRG1 ( INF) 

TAMP ( INF) -AMP ( INF) 

ITRS ( INF) -IRS ( INF) 

2275 CONTINUE 
RETURN 
END 

SUBROUTINE TIMECO ( IDN , TMN, ICNT1 , ICC, IYYDDD) 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
c c 

C TIMECO CHECKS FOR TIME COINCIDENCE BETWEEN DFS TO SEE C 


C HOW MANY DF FIXES SHOULD BE USED IN COMPUTING A SOLN C 



C 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

COMMON/ COM1/ 1 BFLAG , NREC , INCHK , IN (5) , IN2 (5) 

COMMON/ SOLN/SNKA , SSGN L , TLAT , TLON , RADIUS , AREA , TSOLN , 

* I YDSN , I YR , MON , NDAY , NRS , NBR , IRES , NEVT , SMA , SMI , ORIEN 
COMMON/ SUMRY/KTWOER , NTHRS (4 ) , NBADF , NBDF2 , NDUP ( 4 ) , NOVR ( 4 ) , 
*NDFHT(4) , KDVG 

COMMON/ KOUNTR/ KNT 1 , KNT2 , KNT3 , KNT4 , KNTHH , MTKNT , IDFTST , ISECOD 
DOUBLE PRECISION LTIM 
DIMENSION IYYDDD(4) 

TCO-O. 016 
TCP-0.008 
C 

C. . FIRST TIME THROUGH HERE ONLY 
C ICNT1 IS NUMBER OF DFS IN TIME COINCIDENCE 

C 
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IF (INCHK. EQ. X) GOTO 15 
IF (IN(IDN) .BQ.O)THEN 
ICC-0 
IN(IDN)-1 
IKCHK-1 
FTIM-TMN 
LTIM-FTIM 

PTTM*T'TTM+TPP 
CLOCK-FTIM+TCO 
LIDF-IDN 
ICNT-1 
GOTO 20 
ENDIF 

SUBSEQUENT PASSES GO THRU HERE 

CHECK FOR CHANGE OF DAY AFTER MIDNIGHT 

IF(IYYDDD(IDN) . GT. IYYDDD(LIDF) ) THEN 
IF ( (TMN+240000. ) .GT. CLOCK) ITMN-1 
IF ( ITMN . EQ . 1 ) THEN 
TMN-TMN+ 2 4 0 0 0 0 . 

WRITE (6,202) TMN, LTIM 

FORMAT (IX, 'DAY-TIME CHANGE CUR TIME-' , F12 . 3 , 2X, 'LAST-' , 
* F12.3) 

ENDIF 

ENDIF 


C 

C. . SOLN EXISTS WITHIN THE TIME CORRELATION WINDOW 
C 

IF (TMN.GT. CLOCK) THEN 
IF ( ICNT . GE . 2 ) THEN 
IN (5) —1 
ICNT1— ICNT 
TSOLN-FTIM 
IYDSN— IYYDDD(LIDF) 

ENDIF 
GOTO 30 
ENDIF 
C 

c. . CURRENT TIME STILL WITHIN CORRELATION WINDOW 
C 

IF (TMN. LE . CLOCK) THEN 
IF (TMN . GT . CTIM) THEN 
IF ( ICNT . EQ . 1 ) GOTO 30 
C 

C.. CHECK FOR DUPLICATE DF IN TIME WINDOW 
C 

IF (IN (IDN) . EQ . 1 ) THEN 
35 IN (5) -1 

ICNT1— ICNT 
TSOLN-FTIM 
IYDSN-IYYDDD(LIDF) 

GOTO 30 
ENDIF 
ENDIF 


C 

C. . 


DO NOT ALLOW TIME TO DECREASE FROM EVENT I TO EVENT 1+1 


C 

I F ( TMN . LT . LTIM ) THEN 
DO 40 I— 1 , 4 
IN ( I ) —0 

40 CONTINUE 

C WRITE(6, 300) TMN , IDN, LTIM, LI DF 

C300 FORMAT (5X, 'CURRENT TIME ',F12.3,' FOR DF',12,' 

C *' LAST TIME OF ',F12.3,' FOR DF',12) 


LESS THAN', 


C MTKNT IS NUMBER OF OCCURRENCES OF CURRENT TIME LESS THAN PREVIOUS 
C 

MTKNT- MTKNT+1 
GOTO 30 
ENDIF 

C.. LOOK FOR MULTIPLE SOLNS WITHIN TIME WINDOW 

C IF(TMN.LE.CTIM.AND.ICNT.GE.2.AND.IN(IDN).EQ.l) GOTO 35 

IF (TMN. LE. CTIM. AND. ICNT. EQ. 1. AND. IDN, EQ.LIDF) GOTO 30 
C 

ICNT— ICNT+1 
IN (IDN) —1 

C. . RECOVER THE PREVIOUS DF AND PUT IN SOLN SPACE 
C 

I F ( ITMP . EQ . 1) THEN 
ITMP— 0 



poor 


•’ /V ^ $ «r> 

Ot'AL’Tr 


159 



c 

c. . 

c 

30 


C20 

c 

C23 

C 

c 

c 

20 


IN(LIDF)-1 
ENDIF 
LTIM-TMN 
LIDF"IDN 
GOTO 20 
ENDIF 

reset the pointers and tike variables 

ITKP-1 

FTIM-TMN 

LTIM-FTIM 

CTIM-FTIM+TCP 

CLOCK* FTIM+ TOO 

LIDF-IDN 

ICNT-1 

IF (CLOCK. GT. 201815 . .AND. CLOCK. LT. 201818 . )THEN 
WRITE (6,23) ICNT1 , IDN, IN, ICNT, LIDF, LTIM , FT IK, CLOCK, CTIM 
FORMAT ( IX, ' ICNT1, IDN, IN, ICNT, LIDF ' , 9I6/1X, ' LTIM, FTIM, CLOCK. 
* CTIM FOR S23' , 4F12 . 3) ' ' ' 

ENDIF 
RETURN 
RETURN 
END 

SUBROUTINE ERSTAT(0BSA2 ,RXM, DIFDEG, CBRG2) 


C 

C 

c 

c 


SUBROUTINE ERSTAT SORTS DF ANGLE ERROR AS A FUNCTION OF 
ANGLE FOR ALL FIXES, AND IN BINS WITH ERROR ELLIPS LESS 
THAN 10 KM AND FLASHES WITHIN 400 KM OF DF. 


C C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

C 

IMPLICIT DOUBLE PRECISION (A-H, 0-2) 

COMMON/ COM1/ IBFLAG , NREC , INCHX , IN ( 5 ) , IN2 ( 5 ) 

COMMON/ SOLN/ SNKA, SSGNL, TLAT, TLON , RADIUS , AREA, TSOLN, 

*1 YDSN , I YR , MON , NDAY , NRS , NBR , IRES , NEVT , SMA , SKI , ORIEN 
COMMON/ DFSTUF/DFE1 (128) ,NS1 (128) , DFR10 (128) ,NS10 (128) , 
ADFR15 (128) ,NS15(12S) ,SQ1(128) ,SQ10(12S) ,SQ15(128) 

COMMON/ SUMRY/KTWOER , NTHRS ( 4 ) , NBADF , NBDF2 , NDUP ( 4 ) ,NOVR(4) . 
*NDFHT ( 4 ) , KDVG ' *' 

COMMON / KOUNTR/ KNT 1 , KNT2 , KNT3 , KNT4 , KNTHH , MTKNT, IDFTST, ISECOD 
COMMON/ECNST/TMX , RSML1 , RSML2 , RBIG , SBIG , BINS I Z 
DIMENSION AZM(129) , ANWAZM (61), CBRG2 ( 4 ) 

DATA AZM/0. ,2. ,5. ,8. r ll., 14., IS. ,19. ,22. ,25. ,28. ,30. ,33., 
*36. ,39. , 42 45. ,47. , 50. , 53. ,56. , 59. ,61. ,64. ,67. ,70. ,73* , 
*75., 78., 81., 84. ,87. ,90., 92. ,95. ,98. , 101. ,104. , 106. , 109. , 
*112. ,115. ,118. ,120. ,123. ,126., 129. ,132. ,135. ,137. ,140. , 
*143., 146., 149., 151., 154., 157., 160., 163., 165., 168., 171., 
*174., 177., 180., 182., 185., 188., 191., 194., 196., 199., 202., 205., 
*208., 210., 213., 216., 219., 222., 225., 227., 230., 233., 236., 239., 
*241., 244., 247., 250., 253., 255., 258., 261., 264., 267., 270., 272., 
*275. ,278. ,281. ,284. ,286. , 289 . , 292 . , 295. , 298 . , 300. , 303 . , 306. , 
*309. ,312. ,3 15., 317. ,320. ,323. ,326. ,329. ,331. ,334. ,337. ,340. , 
*343. ,345. ,348. ,351. ,354. ,357. ,360./ ' 

DATA ANWAZM/0 .,6. ,12. ,18. ,24. ,30. ,36., 42 . , 48 . , 54. , 60. , 

*66., 72., 78., 84., 90., 96., 102., 108., 114., 120., 126., 132., 138., 

* ' 150 * ' 156 * ' 162 * ' 168 • ' 174 • » 18 °* # 186. , 192. , 198 . , 204 . ,210. , 
*216. ,222. ,228. ,234. ,240. ,246. ,252. ,258. ,264. ,270. ,276. ,282. , 

* 288 • ' 2 **‘ f300. ,306. , 312. , 318. ,324. ,330. ,336. ,342. ,348. ,354. 
*360 + / 

LI2-128 
LI 2 "60 
LIM1-LI2-1 

DO 15 II"1 , LI 2 

IF(OBSAZ .GE.AZM(II) . AND. OBSAZ . LT. AZM(II+1) )THEN 
IF (OBSAZ.GE. ANWAZM (II) . AND. OBSAZ. LT. ANWAZM (1 1+1) ) THEN 
DFS1(II)-DFE1(II)+DIFDEG 
SQ1 (II) -SQ1 (II) +DIFDEG**2 
NS1 (II) »NS1 (II) +1 

IF (ABS (DIFDEG) . LE.TMX) THEN 

IF(SMA.GT.0. . AND. SMA. LE. SBIG. AND. RKK.GE.RSNL1. AND. 

* RKM.LT. RBIG) THEN 

DFR10 (II) "DFR10 (II) +DIFDEG 
SQ10 (II) -SQ10 (II) +DIFDEG**2 
NS 10 (II) "NS 10 (II) +1 
ENDIF 

IF(SMA. GT . 0 . .AND. SMA. LB. SBIG. AND.RKM. GE. RSML2 . AND. 

* RKM.LT. RBIG) THEN 

DFR1 5 (II) "DFR15 (II) +DI FDEG 
SQ15(II)-SQ15(II)+DIFDEG**2 
NS 15 (II) "NS 15 (II) +1 
ENDIF 
ENDIF 
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C IF (II . EQ. 14) THEN 

C WRITE (6,40) II , OBSAZ , DIF DEG, ANWAZM (II) ,DFE1(II) , NS1 (II) , SMA, 

C * DFRIO(II) , NS 10 (II) , DFR15 (II) ,NS15(II) ,SQ1(II) ,SQ10(II) , 

C * SQ15 (II) 

40 FORMAT (IX, 'ERR STATS ' , 14 , 4F8 . 2 , 14 , 2F8 . 2 , 14 , F8 . 2 , 14 , 3 (F8 . 2 , 

* 2X) ) 

WRITE (6,44) TSOLN , TLAT , TLON , IN , IN2 , CBRG2 
FORMAT (IX, 'SOLN ' , F12 . 3 , 2F10 . 4 , 1012 , 4F10 . 4) 

ENDIF 

ENDIF 
CONTINUE 


C 
44 

C 

c 

15 
C 

RETURN 

END 

SUBROUTINE ERSUM 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c c 

c SUMMARIZES AND PRINTS THE DF ERROR STATISTICS C 

c c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

IMPLICIT DOUBLE PRECISION (A-H , O-Z) 

CHARACTER GDATE * 8 

COMMON/ SO LN/SNKA , SSGNL, TLAT , TLON , RADIUS , AREA , TSOLN , 

* I YDSN , IYR , MON , NDAY , NRS , NBR, IRES , NEVT , SMA , SMI , ORIEN 
COMMON/ DFSTUF/DFE 1(128) ,NS1(128) ,DFR10(128) ,NS10(128) , 
ADFR15 (128) ,NS15(128) ,SQ1(128) ,SQ10(128) ,SQ15(128) 

COMMON/ SUMRY/KTWOER, NTHRS (4 ) , NBADF , NBDF2 , NDUP ( 4 ) , NOVR ( 4 ) , 
*NDFHT(4) , KDVG 

COMMON/ KOUNTR/KNT1 , KNT2 , KNT3 , KNT4 , KNTHH , MTKNT , IDFTST , ISECOD 
COMMON/ ECNST/TMX , RSML1 , RSML2 , RBIG , SBIG , BINS I Z 
. EQUIVALENCE (NCNT, NBR) 

DIMENSION AZM(129) , ANWAZM (61) 

DATA AMES , NT1 , NT10 , NT15/-999 . 99,0,0,0/ 

DATA AZM/0. , 2. ,5. ,8. , 11. , 14 . , 16. , 19 . , 22 . , 25 . , 28 . , 30. , 33 . , 
*36. ,39. ,42. ,45. ,47. ,50. ,53. ,56. ,59. ,61. ,64. ,67. ,70. ,73. , 
*75. ,78. ,81. ,84 . ,87. ,90. ,92. , 95 . , 98 . , 101 . , 104 . , 106. ,109. , 


*112 . 

,115. 

,118. 

,120. 

,123. ,126. 

,129., 

*143. 

,146. 

,149. 

,151. 

,154. ,157. 

, 160. , 

*174. 

,177. 

,180. 

,182. 

,185. ,188. 

,191., 

*208. 

,210. 

,213. 

,216. 

,219. ,222. 

,225., 

*241. 

,244. 

,247. 

,250. 

,253. ,255. 

,258., 

*275. 

,278. 

,281. 

,284. 

,286. ,289. 

,292., 

*309. 

,312. 

,315. 

,317. 

,320. ,323. 

,326. , 

*343. 

,345. 

,348. 

,351. 

,354. ,357. 

,360./ 


3000 


510 


190 


C 

c. . 
c 


210 


DATA ANWAZM/ 0 . , 6 . , 12 . , 18 . , 24 . , 30. , 36. , 42 . , 48 . , 54 . , 60. , 

*66. ,72. ,78. ,84. ,90. ,96. , 102. ,108. ,114. , 120. ,126. ,132. ,138. , 

*144. , 150. , 156. , 162 . ,168. , 174 . , 180 . , 186 . , 192 . , 198 . ,204. ,210. , 

*216. ,222. ,228. ,234. ,240. ,246. ,252. ,258. ,264. ,270. ,276. ,282., 

*288. , 294. , 300. ,306. , 312 . , 3 18 . , 324 . , 330 . , 336 . , 342. ,348. ,354. , 
*360./ 

LI3-128 
LI 3 “60 

WRITE (6, 3000) NEVT , TSOLN , I YR , MON , NDAY , TLAT , TLON , SNKA , SSGNL, 

* NRS , RADIUS , AREA 

FORMAT (/IX, 'LAST SOLN ON TAPE-' , 16 , 2X, F12 . 3 , IX, 312 , 2F12 . 6 , 2X, 

* 2F10. 2 , 12 , 2F8 . 2/) 

WRITE ( 6 , 510) TSOLN, KNTHH, IYR, MON, NDAY 

FORMAT (//IX, 'FINAL HOURLY SUMMARY ON TAPE ENDING AT ',F12.3, 

*' IS' , 16, 4X, 'DATE-' , 312//) 

WRITE ( 6 , 190) IDFTST, TMX , RSML1 , RSML2 , RBIG, SBIG, BINSIZ 
FORMAT (1H1/ IX, 'ERROR SUMMARY FOR DF',I2//4X,' CONSTRAINTS:', 

*5X, 'MIN. ANGLE DEV. -' , F4 . 1 , ' RANGES-' , 3F6 . 1 , ' MIN. SMA-',F6.1, 

*' BEARING ANGLE BIN-',F4.1/) 

IF IDFTST— 0 THEN NOT WANTING SITE ERRORS, SKIP TO SUMMARY 
IF (IDFTST. EQ.O) GOTO 600 
WRITE (6,210) 

FORMAT (//4X, 'I',4X, 'AZM',4X, 'N',5X, 

*'AV1',6X, 'VAR',8X, ' SD' , 6X, 'N' , 6X, ' AV10' , 6X, 'VAR10' , 6X, 'SD10',6X, 
*'N',6X, ' AV15 ' , 5X, ' VAR15 ' , 6X, 'SD15'/) 

DO 15 IK-1, LI3 

IF(NS1 (IK) .LE. 1)THEN 
Xl-DFEl(IK) 

VI— AMES 
SD1-AMES 
GOTO 300 
ENDIF 

AS1— FLOAT (NS1 ( IK) ) 

ASM 1 -AS 1- 1 . 

XI— DFE1 ( IK) /AS1 



Vl» (SQ1 (IK) - DFE1 (IK) **2/ASl) /ASM1 
SDl-SQRT(Vl) 

300 IF(NS10(IK) .LE.1)THEN 

XIO-DFRIO(IK) 

VI 0" AMES 
SD10-AMES 
GOTO 400 
EKDIF 

AS10-FLOAT(NS10(IK) ) 

ASM10— AS10-1 • 

X10-DFR10(IK)/AS10 

V10-(SQ10(IK) - DFRIO(IK) **2/AS10)/ASM10 
SDIO-SQRT(VIO) 

400 IF (NS15 (IK) . LE. 1)THEN 

X15-DFR15(IK) 

V15-AMES 
SD15-AMES 
GOTO 500 
EKDIF 

AS15-FL0AT(NS15(IK) ) 

ASM15-AS15-1. 

X15-DFR15 ( IK) /AS 15 

V15-(SQ15(IK) - DFR15(IK) **2/AS15)/ASM15 
SD15-SQRT(V15) 

C 

C. . WRITE SITE ERROR INFO TO DISK FILE AND TO PRINTER 
C 

LUIN-10+IDFTST 

500 WRITE ( LUIN ,110) IK, ANWAZM(IK) ,NS10(IK) ,XX0,SDX0 

110 FORMAT ( 14 , IX, F5 . 1 , 15 , 2F10 . 2 ) 

WRITE(6, 100) IK, ANWAZM(IK) ,NS1(IK) ,X1,V1,SD1,NS10(IK) , X10, 

* V10, SD10, NS 15 (IK) , X15, V15 , SD15 

100 FORMAT (1X,I4,2X,F5.1,I5,3F10.2,I5,3F10.2,I5, 3F10. 2) 

NT1-NT1+NS1(IK) 

NT10-NT10+NS10(IK) 

NT15-NT15+NS15(IK) 

15 CONTINUE 

600 WRITE(6, 601) NT1, NT10, NT15 

601 FORMAT (IX, 128 ( ' . ' )//lX, 'TOTAL NO. OF INPUTS 4X, 'Nl-', 16, 4X, 
*'N10-',I6,4X, 'N15-',I6) 

WRITE (6, 611) KNT1 , KNT2 , KNT3 , KNT4 , KTWOER , NBADF , NBDF2 , 

*NDFHT , NTHRS , NDUP , MTKNT , NOVR , KDVG 
611 FORMAT(/2X,'NO. OF 1,2,3, OR 4 DF SOLNS-' , 4 (16, 2X) / 

*2X,'NO. OF BAD REPEAT SOLNS-' ,15,/ 

*2X, 'NO. OF NO-FIXES (NO SOLN) FROM FFIX- ',16,/ 

*2X,'NO. OF BAD DF AZMS, ANGLE DEVIATION FROM SOLN TOO BIG-' 16 / 
*2X, 'TOTAL NO. OF DF HITS-' , 4 (16, 2X)/ 

*2X, 'NO. OF DF HITS BELOW MIN. AMPLITUDE OF 10-', 416,/ 

*2X, 'NO. OF DF DUPLICATE HITS WITHIN TIME CORR. WINDOW-' 416,/ 
*2X,'NO. OF OCCURRENCES WITH LATEST TIME LESS THAN PREVIOUS-' 16/ 
*2X, 'NO. OF DF OVERANGES-' ,416/ 

*2X, 'NO. OF DIVERGENT DF ANGLE PAIRS-', 16/) 

CALL DATEG(GDATE) 

CALL TIME (ITIME) 

IHTIM-INT ( ITIME/ 3 60000) 

AHTIM— FLOAT ( ITIME)/360000 . 

AKTIM— (AHTIM-IHTIM) *60. 

IMTIM— INT ( AMTIM) 

ASTIM-60 . * (AMTIM- IMTIM) 

STIM-FLOAT( (IHTIM*100 + IMTIM) *100) +ASTIM 
WRITE (6, 3) GDATE , STIM, IDFTST 

3 FORMAT(//20X,'END LLP QUALITY CONTROL ANALYSIS '/2 OX, ' ( ' , A8, 2X, 
*F10.2 , ' ) ' , 10X, 'DF TEST— ',12) 1 ' ' ' 

STOP 
END 

FUNCTION FKA(TAMP) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c c 

C FUNCTION FKA RECEIVES THE AMPLITUDES (AMP) OF THE DF'S C 

C AND USES THE GREATEST SIGNAL STRENGTH (BMAX) TO COMPUTE C 

C A RANGE NORMALIZED SIGNAL STRENGTH - DI ST* BMAX, WHERE C 

C DIST IS THE DISTANCE BETWEEN THE DF AND THE SOLN POINT. C 

C THIS PRODUCT IS DIVIDED BY 298, A CALIBRATION FACTOR C 

C BASED ON THE LLP ANTENNA AND UN ET AL'S TRANSMISSION C 

C LINE MODEL OF THE RETURN STROKE, RETURNS PEAK CURRENT C 

C ESTIMATE FOR DISCHARGE C 

C C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

COMMON/ COM1/I BFLAG , NREC , INCHK, IN (5) , IN2 (5) 

COMMON/ SOIM/SNKA , SSGNL, TLAT, TLON , RADIUS , AREA, TSOLN , 

*IYDSN, IYR, MON , NDAY , NRS , NBR , IRES , NEVT, SMA, SMI , ORIEN 
COMMON/DFSTUF/DFE1 ( 128 ) ,NS1(128) , DFR10 (128) ,NS10(128) , 

ADFR15 ( 128) ,NS15(128) ,SQ1(128) ,SQ10 (128) , SQ15 ( 128) 
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c 

c.. 

c 


COMMON/SUMRY/KTWOER,NTHRS(4) , NBADF , NBDF2 , NDUP ( 4 ) ,N0VR(4) , 
*NDFHT ( 4 ) KDVG 

COMMON/ K0UNTR/KNT1 , KNT2 , KNT3 , KNT4 , KNTHH , MTKNT , IDFTST, ISECOD 
EQUIVALENCE (NCNT, NBR) 

DIMENSION TAMP (4) 

DIMENSION ALATI ( 4 ) , AIDNGI ( 4 ) 

DATA ALATI/34. 649167,35. 399167,35. 83750,34. 716667/ 

DATA ALONGI/86 . 669167 , 86 . 076944 , 87 . 443889 , 87 . 881667/ 

DATA DTR,ERAD2,NST/0. 01745329, 6371. ,4/ 

ISPOL“0 

BMAX-0. 

FIND ID AND AMP OF LARGEST SIGNAL STRENGTH 


DO 20 1-1, NCNT 

IF (TAMP ( IN2 ( I ) ) . GT . 0 . ) ISPOL-ISPOL* 1 
IF(ABS(TAMP(IN2(I) )) .CT. BMAX) THEN 
BMAX-ABS (TAMP (IN2 (I) } ) 

KDF«IN2(I) 

ENDIF 
20 CONTINUE 

IF ( KDF . EQ . 0) THEN 

WRITE (6,50) TSOLN , KDF , IN2 , AZ2 ,R2 , BMAX , SPOL, FKA , NCNT, TLAT , TLON , 

* RADIUS , TAMP 
FKA— 9999. 

RETURN 

ENDIF 

C 

C.. FIND POLARITY OF MAX AMPLITUDE (POSITIVE IFF ALL DFS POSITIVE) 

C 

SPOL— 1. 

IF(ISPOL.EQ.NCNT) SPOL-1. 

RMAX2- BMAX* SPOL 
SSGNL-RMAX2 
C 

C.. FIND RANGE NORMALIZED DISTANCE 
C 

SLAT - ALATI (KDF) *DTR 
SLONG - ALONG I (KDF) *DTR 
XTL2- TLAT*DTR 
XTL02- TLON*DTR 

CALL AZRN (SLAT, SLONG, XTL2 , XTL02 , AZ2 , R2 ) 

C 

FKA - R2*ERAD2*RMAX2/298. 

C 

C WRITE(6,50) TSOLN, KDF, IN2,A22,R2, BMAX, SPOL, FKA, NCNT, TLAT, TLON, 

C *RADIUS , TAMP 

50 FORMAT (/IX, 'KA SOU*.. ' , F12 . 3 , IX, 614 , 5F12 . 4/1X, ' NCNT' , 12 , 7F10. 4) 
RETURN 
END 

FUNCTION SITERR(IDF, BRG) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


C c 
C FUNCTION SITERR RECEIVES A DIRECTION FINDER ID AND C 
C AN UNCORRECTED BEARING HAVING A KNOWN SITE ERROR WHICH C 
C CAN BE FOUND BY A LINEAR INTERPOLATION PROCEDURE APPLIED C 
C TO A LOOK-UP TABLE FOR THAT DF C 
C C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

COMMON/ DFSTUF/DFE1 (128) ,NS1(128) ,DFR10(126) ,NS10(128) , 
*DFR15 ( 128) , NS15 ( 128) ,SQ1(128) ,SQ10(128) ,SQ15(128) 

COMMON/ KOUNTR/KNT1 , KNT2 , KNT3 , KNT4 , KNTHH, MTKNT, IDFTST, ISECOD 
DIMENSION DF1 (60) , DF2 (60) , DF3 (60) ,DF4(60) ,AZM(61) 

DIMENSION DFA(128) ,DFB(128) ,DFC(128) ,DFD(128) , AZMI (129) 

DATA DF1/-2 .47,-1.89,-1.65,-0.13,-1.02,-2.69,-3.57, 

*-3.45, -4. 00, -2. 78, -7.0 ,0.83,1.08,0.89, 

*2.14,2.69,3.4,4.09,4.90,5.58,5.39,5.38,5.59,5.89,6.03,6.72, 
*6.51,6.44,6.86,7.14,7.04,6.56,6.29,6.21,-4.2,5.13, 
*4.61,4.85,5.05,3.50,3.46,2.67,2.72,1.68,-1.2,-1.21,-1.20, 
*-0.6, -0.42, -0.3, -0.4, -.87, -1.35, -1.3, -1.65, -4. 18, 

*-4.5, -3. 99, -4. 29, -3. 53/ 

DATA DF2/-2. 1,-1. 4, -1.6, -1.6, -1.1, -0.2, 0.01, -.15, -.21, 

*-.46, -.74, -.61, .96, 0.2, -1.7, -3. 8, -5. 5, -6. 1,-6. 6, -7. 2, 

*-5.6, -5. 0,-6. 8, -9. 9, -12. 2, -14. 0,-14. 8, -13. 5, -11. 6, -9. 9, 
*—8.0,— 6.5,— 4.6,— 3.2,— 0.8,— 3.9,2.8,2.50, 1.8,1. ,—6.4, 1.5, 

*2 . 50 , 3 . 9 , 5 . 4 , 6 . 8 , 3 . 00 , -2 . 1 , -3 . 2 , -2 . 8 , -3 . 7 , -6 . 6 , -9 . 70 , -10 . 7 , 
*-8.5, -6. 1,-5. 3, -3. 8, -2. 8, -2. 3/ 
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DATA DF3/-1. 5, -1.2, -2. 0, -2. 0,-0. 6, -1.7, -1.4, -1.3, -1.2, 

#-l. 2,-1. 1,-0. 9,-0. 3, . 1, 0.3,0.50,-. 30, -1.9, -0.8, 0.4, 

♦•0.2 ,-1.5,— 2. 4, -3. 3, -4. 3, -5 .2,-6 .1,-7. 0,-8. 0,-8. 9, -9 . 8 , 

-8.9, -8 .0,-7. 1,-6. 2, -5. 3, -4. 4, -5. 5, -2. 0,-1. 5, -1.0, -0.46, 

0.05, 0.6, 1.1, 1.6,2. 1,2. 6, 3. 1,3. 1,2. 9, 2. 1,0. 2, -.7, -0.8, 

-2. 7, -2. 8, -1.5, -1.3, -1.4/ 

DATA DF4/— • 4,— .7,— .6, .24, . 67 ,1.9, 2.0, 2.0, 2.2, 2*3, 2. 2, 

1.8. 5. 2. 4. 2. 4. 1.4. 7. 5. 1.5. 2. 5. 6. 6. 4. 6. 4. 5. 9. 5. 0. 4,7, 4. 4, 4.1, 

3. 7. 3. 2. 3. 0. 2. 7. 1.9. 1.0. .22.. 64. 1.3. 2. 5. 2. 6.2. 3. 1.9. 1.3, 

.68, .56, .62, .15, .07, -.2, .07, .05, .10, -.2, -.6, -.3, .1, .5,1.3, 

1.7. 1.2. . 1,-0. 01, -.3/ 

DATA DFA/-3 .7, -2. 3, -1.5, -2.0, -.9, . 5, . 9 , . 4 , -.7 , -1 . 3 , 

—1.4 , —1.9, —1 . 8, —3 • 6, —4 .4, —5. 5,— 4 .6,-3 . 6, —4 .2, —2. 8 , 

—4 .2,— 2.7, — *2,— 1.7,— 1.1, 1*4, —2.2, 1.0, .5, .9, 

•1,2., .8,0.0, *4, 1.5, 1*5, 3. 3, 1.7, 3.2, 

2. 7. 1.8. 2. 9. 2. 5. 2. 8. 3. 0. 3. 3. 4. 2. 5. 3. 4. 8, 

4. 9. 3. 9. 4. 3. 3. 5. 4. 5. 4. 4. 5. 1.6. 0. 6. 4. 6. 8, 

7. 9. 8. 9. 8. 9. 9. 4. 6. 1.6. 0. 5. 3. 5. 6. 4. 7. 4. 8, 

3. 8, 3. 4, 4. 3, 3.9, 4. 1,4.6, 5. 6, 5. 3, 5. 6, 6.1, 

4. 7. 4. 3. 5. 8. 6. 0. 3. 6. 3. 6. 3. 7. 2. 0.2. 8, -3. 4, 

2. 0,-1. 5,. 5, -4.0,. 4, -3. 7, -2. 2, -3. 8, -1.2, -4. 4, 

-3. 0,-2. 6, -4. 2, -1.7, -5. 2, -4. 9, -4. 5, -7. 1,-11. 1,-5. 4, 

-7. 2, -5. 8, -2. 8, 1.2, 1.5, -2. 4, -1.3, -1.5, -3. 8, -2. 5, 

—5 • 4 , — • 2 , —2 . 0 , — . 8 , —1 .2,— .8,— .4,1. 5/ 

DATA DFB/-2 .3, -2. 3, -1.6, -1.3, -1.5, -1.6, -1.7, -1.7, -1.5, 

1.2, -0.8, -0.5, -0.3, -0.1,0. ,0. ,-0.1, -0.2, -0.2, -0.3, 

0.4, -0.6, -0.7, -0.8, -0.7, -0.4, 0.1, 0.7, 1.0, 0.5, -0.6, -1.5, 

-2. 3, -3. 3, -4. 5, -5. 3, -5. 7, -6. 0,-6. 2, -6. 5, -6. 9, -7. 2, -7.0, 

6. 5, -5. 8, -5. 2, -5. 0,-5. 3, -6. 2, -7. 7, -9. 4, -10. 8, -11. 9, -12. 9, 
13.8,-14.4,-14.7,-14.9,-14.7,-13.9,-12.8,-11.9,-11.1, 

-10 .2 , —9 .2,— 8*3,— 7.5,— 6.8,— 5.9,— 5.0,— 4.3,— 3.5,— 2.4,— 1.2, 

0.,1. ,1.9, 2. 5, 2. 8, 2. 8, 2. 6, 2. 3, 1.9, 1.5, 1.1, 0.9, 0.8,1., 

1.4. 1.7. 2. . 2.3.3. .3. 7. 4. 3. 5. 1.6. 1.6. 8. 6. 1.3. 9. 1.1, -1.2, 

-2. 8, -3 .3,-3. 1,-2 .9,-2 .7, -2. 8, -3 .4, -4. 6, -6. 1,-7. 6, -9.1, 

-10.2, -10 . 7 , -10 . 3 , -9 . 0 , -7 . 3 , -6 . 3 , -5 . 8 , -5 . 4 , -5 . , -4 . 6 , -4 . 1 , 

-3. 5, -3. ,-2.3, -2. 3/ 

DATA DFC/-1. 4, -1.4, -1.6, -1.3, -1.1, -1.6, -2. 4, -2. 3, -1.4, 

— 0 .7,— 0.6,— 1.2,— 1.7,— 1.6,— 1.4,— 1.3,— 1.3,— 1.3,— 1.2,— 1.1, 

— 1.2,— 1. 2 , — 1 . 1 , -1 • 0, — . 9, — • 8 , — .7 , — . 4 , — . 1, .1, .2, .3, .4, .5, 

.5,0. ,-.9, -1.7, -1.9, -1.2,-. 1, .4, .4,-. 1,-. 9, -1.6, -2.1, 

-2. 5, -2. 7, -2. 9, -2. 9, -2. 9, -2. 8, -2. 8, -2. 8, -2. 9, -3. 1,-3. 4, 

-3 . 7 , -4 . , — 4 . 1 , — 4 . 2 , — 4 . 3 , — 4 • 4 ,— 4 .5, -4 . 5 , -4 . 5,— 4. 5, — 4 . 5, 

-4. 5, -4. 1,-3. 6, -3. 5, -4. 1,-4. 5, -4. 4, -4. ,-3.4, -2. 7, -2.1, 

—1 *7,— 1.4,— 1.1,— 1.1,— 1.2,— 1.3,— 1.4,— 1.5,— 1.3,— .9,-1., 

-1.5, -1.9, -2. 1,-2. ,-1.7, -1.4, -1.1, —.6, .3, 1.4, 1.6, 1.4, 

1.. 1.3. 2. 6. 3. 8. 3. 9. 3. 2. 2. 6. 2. 2. 1.6. .6, -.2, -.6, -.8, -.8, 

-1.3, -2. 4, -3. 1,-2. 9, -2. 5, -2. 3, -1.8, -1.2, -1.2, -1.4, -1.4/ 

DATA DFD/-.3,-.3, -.6, -.7, -.7, -.7, -.4, . 1, . 5 , . 6, . 9, 1 . 4 , 

1.8. 2.. 2.. 2.. 2.. 2. 1.2. 2. 2. 3. 2. 3. 2. 3.2. 2. 2. 1.1. 9. 1.6. 1.7, 

2. 6, 3. 6, 4. 1,4. 1,4. 1,4. 3, 4. 6, 4. 9, 5. 1,5. 2, 5. 2, 5. 3, 5. 5, 5. 9, 

6. 3. 6.5. 6. 6. 6. 5. 6. 3. 6. . 5. 6. 5. 2. 4. 9. 4. 7. 4. 6. 4.5. 4.3. 4.2. 4., 

3.8. 3. 6. 3. 5. 3. 3. 3.1. 3. . 2. 9. 2. 7. 2. 5. 2.1. 1,7, 1.2, .7, .3, .2, .5, 
.9, 1.2, 1.6, 2. 1,2. 4, 2.6, 2.6,2. 5, 2. 4, 2.2, 2., 1.7, 1.4, 1.1, .8, .5, 
.5, .7, .8, .7, .4, .2, .1, • 1, • • 1, * .2, -.1,0. ,.2, .l,0.,-.l,— .1,— .2, 

— .3,— .4,— .6,~*6,— .4,0. , • 1 , • 1 , .4, .8, 1.2, 1.5, 1.7, 1.6, 1.3, .9, 

.5, .2, .2 , .1,— .3,— . 3/ 

DATA AZM/0. ,6. , 12. , 18. , 24 . , 30. , 36 . , 42 . , 48 . , 54 . , 60. , 66. ,72 . , 

78 .. 84. .90. .96. .102. .108.. 114 , ,120. ,126. ,132. ,138. ,144. , 150. , 

156.. 162. .168. .174.. 180.. 186.. 192. .198.. 204.. 210.. 216.. 222., 

228. . 234. .240. .246. .252. .258. .264. .270. .276. .282.. 2 88., 

294. . 300. .306. .3 12 . . 3 18 . .324. .330. .336. .342. . 348 , , 354 • , 360 •/ 
DATA AZH1/0. ,2. ,5. ,8. ,11. ,14. ,16. ,19. ,22., 25. , 26., 30., 33. , 

36 . . 39 ..42. .45. .47.. 50. . 53 . . 56 , , 59. , 61 .,64. ,67. ,70., 73 • , 

75.. 78.. 81.. 84.. 87.. 90.. 92.. 95.. 98.. 101.. 104.. 106.. 109., 

112. . 115. .118. .120. .123. .126.. 129. .132. .135. .137. .140. , 

143 . . 146 . . 149. . 151 . . 154 . . 157. . 160. . 163 . . 165. . 168. . 171. , 

174.. 177.. 180.. 182.. 185.. 188.. 191.. 194.. 196.. 199.. 202.. 205., 

208.. 210.. 213.. 2 16.. 2 19.. 222.. 225.. 227.. 230.. 233.. 236.. 239., 

241.. 244. .247. .250. .253.. 255. .258.. 261.. 264.. 267. .270. .272. , 

275.. 278.. 281. .284. .286.. 289.. 292.. 295.. 298.. 300.. 303.. 306., 

309.. 312.. 315. .317. .320.. 323. .326. .329. . 331. .334. .337. .340. , 

343.. 345.. 348. .351.. 354.. 357.. 360,/ 

DTR-0. 01745329 

LI-128 
LI-60 
HI— 1*1-1 

I F ( IS ECOD . EQ . 3 ) THEN 
WRITE (6, 32) 

FORMAT (1H1/ IX, 'SITE ERROR TABLE'/) 

DO 5 I— 1,LI 

WRITE (6,2) I,AZM(I) ,DF1(I) , DFZ(I) ,DF3(I) ,DF4(I) 

FORMAT ( IX , 14 , 5F10 . 2) 

CONTINUE 
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ISECOD— 1 
WRITE (6, 33 ) 

33 FORMAT (1H1) 

ENDIF 

NDF-4 

DO 10 I*»l , LI 

I F ( BRG . GE . AZM ( I ) . AND. BRG. LT. AZM (1+1) )THEN 
C 

IF ( IDF . EQ . 1 ) THEN 
ERl-DFl(I) 

ER2*DF1(I+1) 

AZl-AZM(I) 

AZ2-AZM(I+1) 

ELSE IF ( IDF . EQ . 2 ) THEN 
ER1«DF2(I) 

ER2«DF2(I+1) 

AZ1=AZM (I) 

AZ2-AZM(I+1) 

ELSE I F ( I DF . EQ . 3 ) THEN 
ER1-DF3(I) 

ER2«DF3(I+1) 

AZl-AZM(I) 

AZ2"AZM (1+1) 

ELSE IF (IDF. EQ. 4) THEN 
ER1-DF4(I) 

ER2-DF4 (1+1) 

AZl-AZM(I) 

AZ2-AZM(I+1) 

ENDIF 

ENDIF 

10 CONTINUE 

C 

C. . PERFORM INTERPOLATION 

C 


SLOPE" (ER2-ER1) / (AZ2-AZ1) 
SITERR"ER2 -SLOPE* (AZ2-BRG) 

C IF(IDF.GT.O) GOTO 909 

C 

C. . TEMPORARY CORRECTION TO DFS3 AND 1 


C 

C IF(IDF . EQ. 3 ) THEN , 

C SITERR— 3 . 5+ . 6*SIN ( 2 . *BRG*DTR) -5 . l*COS (2 . *BRG*DTR) 

C ENDIF 

IF(IDF.EQ.IDFTST) SITERR-0. 

C IF (IDF. EQ. 1) SITERR-0. 

C WRITE (6, 25) IDF, BRG, SITERR 

25 FORMAT ( IX, 12 , 2X, 2 (F6 . 2 , 2X) ) 

C DO 120 I-1,LI 

C WRITE (6,35) I,DF2(I) , DF3 (I) , DF4 (I) ,AZM(I) 

35 FORMAT (IX, 14, 4 (F6. 2, 2X) ) 

C120 CONTINUE 

C WRITE (6 , 399) 

399 FORMAT ( IX , 9 GOING BACK TO DFDATA') 

909 RETURN 
END 

s^ccccccccccecc^ cccccccccccccccccccxccccccCTrecccccccccc 


c FUNCTION SITERC IS A MOD TO SITERR USING CONSTRAINT C 
C POINTS FROM VISUAL AND RADR CONFIRMATION OF CELL £ 
C LOCATION. USES A LINEAR INTERPOLATION PROCEDURE APPLIED C 
C TO A LOOK-UP TABLE FOR THAT DF £ 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


c 


■OMMON /SOLN/SNKA , SSGNL , TLAT , TU3N , RADIUS , AREA , TSOLN , 

YDSN , I YR , MON , NDAY , NRS , NBR , IRES , NEVT , SMA , SMI,ORIEN 
COMMON/ DFSTUF/DFE1 ( 128 ) ,NS1(128) ,DFR10(128) ,NS10(128) , 

>FR15 ( 128) ,NS15(128) ,SQ1(128) ,SQ10(128) TO _ nrin 

:OMMON / KOUNTR/ KNT 1 , KNT2 , KNT3 , KNT4 , KNTHH , MTKNT , I DFTST , ISECOD 
)IMENSION DF1 (60) , DF2 (60) , DF3 (60) , DF4 (60) , AZM (61) 

)ATA DF1/4.0, 4 .0, 1 . 65 , -0 . 13 , -2 . 0 , -5 . 0 , 2 . 0 , 
l. 0,-4. 00,-5. 0,-5. 5 ,0.83,1.08,0.89, 

>.14,2.69,3.4,4.09,4.90,5.58,5.39,5.38,5.59,5.89,6.03,6.72, 

i! 51, 6. 44, 6. 86, 7. 14, 7. 04, 6. 56, 6. 29, 6. 21, -4. 2, 5. 13, 

1.61, 4. 85, 5. 05, -1.0, -1.0, 2. 67, 2. 72, 1.68, —1.2, “1.21, 1.20, 

•0.6, -0.4 2, -3.0, -5.0, -.87, -1.3 5, -1.3, -1.2, “.5, 

3ATA DF2/-2*. 1,-1. 4, -1.6, -1.6, -1.1, -0.2, 0.01, -.15, -.21, 

-.46, -.74, -.61,-1. ,“1.5, -1.7, -3. 8, -5. 5, -6. 1,-6. 6, -7. 2, 

-5. 6. -5. 0,-6. 8, -9. 9, -12. 2, -14. 0,-14. 8, -13. 5, -11. 6, -9. 9, 


00147223 

00147323 

00147423 

00147523 

00147723 

00147823 

00148023 
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UUOUM uuouuuuu 


*-8. 0,-6. 5, -4. 6 , -3. 2, -0.8, 0.0, 0.0, 0.00, -2. , -2., -1.5, -4. , 

*-10. , 3. 9, 5. 4, 6. 8, 3. 00 r -2. 1,-3. 2,-2. 8, -3. 7, -6.6, -9. 70, -10. 7, 
*-8.5, -6. 1,-5. 3, -3. 8, -2. 8, -2. 3/ 

DATA DF3/-1. 5, -1.2, -2. 0,-2. 0,-0. 6, -1.7, -1.4, -1.3, -1.2, 

•-1. 2, -1.1, -0.9, -0.3, .1,0. 3, -1.0, -2. 4,-1. 9,-0. 8,-1. , 

*— 3 . 0,-2. 5, 0.0, 0.0, -2.0, -8.0, -12 . , -11 . , -10. , -10. , -9 . 8 , 

*—5.9 ,-5. 0,-4. 5, -4. 0,-3. 5, -3. 0,-2. 5 ,-2. 0,-1. 5, —1.0, —0.46, 
*0.05, 0.6, 1.1, 1.6, 2. 1,2. 6, 3. 1,3. 1,2. 9, 2. 1,0. 2, -.7, -0.8, 

*-2.7, -2. 8, -1.5, -1.3, -1.4/ 

DATA DF4/-.4,-.7,-.6, .24,. 67, 1.9, 2. 0,2. 0,2. 2, 2. 3, 4.0, 

*2. 0,1. 5, 3. 2, 4. 1,4. 7, 5. 1,5. 2, 5. 6, 6. 4, 6. 4, 5. 9,4. 5, 2. 0,4. 4, 4.1, 
*3. 7, 3. 2, 3. 0,2. 7, 1.9, 1.0, . 22,. 64, 1.3, 2. 5, 2. 6, 2. 3, 1.9, 1.3, 

*.68, .56, .62, .15, .07, -.2, .07, .05, .10, -.2, -.6, -.3, .1, .5,1.3, 
*1.7, 1.2, .1,-0. 01, -.3/ 

DATA AZH/0. , 6 . , 12 . , 18 . ,24 . ,30. , 36. , 42 . , 48 . ,54 . , 60. , 66. , 72 . , 
*78. ,84. ,90. ,96. ,102. ,108. ,114. ,120. ,126. ,132. ,138., 144. ,150. , 
*156. ,162. ,168. ,174. ,180. , 186 . , 192 . , 198 . , 204 . , 210 . , 216. ,222. , 
*2 28. ,23 4. ,240. ,246. , 252 . , 258 . , 264 . ,270. ,276. ,282., 288., 

*294., 300. ,306. ,312. ,318. ,324., 330. , 336. , 342. , 348. ,354 . , 360 ./ 

I F ( IDF . EQ . IDFTST ) THEN 
SITERC-0 . 0 
RETURN 
END IF 
SLOPE— 0 . 

DTR-0. 01745329 

LI-60 

LX1-LX-1 

IP ( ISECOD . £Q . 1 ) THEN 
WRITE (6,32) 

32 FORMAT(lHl/lX, 'SITE ERROR TABLE'/) 

DO 5 1—1,11 

WRITE(«,2) X,AZH(X) ,DP1(I) ,DF2(I) ,DF3(I) ,DP4(Z) 

2 FORMAT(lX, I4,5F10.2) 

5 CONTINUE 

ISECOD— 1 
WRITE (6,33) 

33 FORMAT (1H1) 

ENDIF 

NDF-4 

DO 10 X-1,LI 

IFfBRG.GE. AZM(I) . AND. BRG.LT.AZWfI+1) )THEN 
C 

IF ( IDF . FQ . 1 ) THEN 
ER1— OF 1(1) 

ER2— OF1 (1*1) 

A21—A2M (I) 

AZ2— AZM (1*1) 

GOTO 20 

ELSE IF ( IDF . EQ . 2 ) THEN 
ER1— 0F2 ( I ) 

ER2-DF2 ( 1+1) 

A21-AZMCI) 

AZ2-AZM(I+1) 

GOTO 20 

ELSE IF ( IDF . EQ . 3 ) THEN 
ER1— OF3{I) 

ER2— 0F3 ( 1+1) 

AZ1— AZM(I) 

AZ2— AZM (1*1) 

GOTO 20 

ELSE IF ( IDF . EQ . 4 ) THEN 
ER1-DF4 (I) 

ER2-DF4 (1*1) 

A21-AZM (I) 

AZ2— AZM (1*1) 

GOTO 20 

ENDIF 

ENDIF 

10 CONTINUE 

PERFORM INTERPOLATION 


0 DELA— AZ2-AZ1 

SLOPE- ( ER2 -ER1 ) /DELA 
SITEROER2 -SLOPE* (A22-BRG) 

IF (IDF . GT. 0) GOTO 919 

TEMPORARy CORRECTION TO DFS3 AND 1 


ORIGINAL PAGE IS 
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IF ( IDF . EQ. 3 ) THEN 

IF (BRC.CZ. 150. . AND. BRG.LE. 180.) THEN 
5XTERC— SITERC-3 
ENDIF 
ENDIF 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

25 

c 

c 

35 

C120 

C 

399 

919 


IF (IDF. EQ. 2) THEN 

IF (BRG. GE. 230 . . AND. BRG. LE . 250 .) THEN 
SITEROSITERC-6 . 

ENDIF 

ENDIF 

IF (IDF . EQ . 1)THEN 

IF (BRG. GE. 330. .AND. BRG, LE. 340. ) THEN 
SITERC-SITERC+3 . 

ENDIF 

ENDIF 

IF ( IDF . EQ . 4 ) THEN 

IF ( BRG. GE. 65 . . AND. BRG . LE. 75 . ) THEN 
SITERC“SITERC+1 . 5 
ENDIF 
ENDIF 

IF ( IDF . EQ . 3 ) THEN 

SITERR-~3 . 5+ . 6*SIN ( 2 . *BRG*DTR) -5 . 1*C0S ( 2 . *BRG*DTR) 
ENDIF 

WRITE (6,25) IDF , BRG, SITERR 
FORMAT (IX, 12, 2X, 2 (F6.2,2X) ) 

DO 120 1=1 , LI , % 

WRITE (6,35) I , DF2 ( I ) ,DF3(I) ,DF4(I) ,AZM(I) 

FORMAT (IX, 14,4 ( F6 . 2 , 2 X ) ) 

CONTINUE 
WRITE(6, 399) 

FORMAT ( IX, ' GOING BACK TO DFDATA' ) 

RETURN 

END 

FUNCTION SITER2 (IDF , BRG) 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c c 

c FUNCTION SITER2 RECEIVES A DIRECTION FINDER ID AND C 

C AN UNCORRECTED BEARING HAVING A KNOWN SITE ERROR WHICH C 

C CAN BE FOUND BY FITTING A POLYNOMIAL OF DEGREE 6. C 

C C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C 

C0MM0N/DFSTUF/DFE1 ( 128 ) ,NS1(128) ,DFR10(128) ,NS10(128) , 

*DFR15 ( 128 ) , NS15 ( 128 ) ,SQ1(128) ,SQ10(128) ,SQ15(128) 

COMMON/ KOUNTR/KNT1 , KNT2 , KNT3 , KNT4 , KNTHH, MTKNT, IDFTST, ISECOD 
DIMENSION DFA ( 60) ,DFB(60) ,DFC(60) ,DFD(60) ,AZM(61) 

DIMENSION AO (4) , ADI (12) ,AD2 (12) , AD3 (12) , AD4 (12) 

DATA AO/ 1.63, -3. 33, -1.47, 1.94/ 

DATA AD1/0 .8788,-4 . 9730,-0.4359, 0. 2113, -0.2401, 0. 2218, 
*0.4777, 0.0991,0.8518, 0. 1006, 0.4715,0. 4008/ 

DATA AD2/-2. 1636, 1.6295, 4. 9684, -2. 5906, -0.0242, 1.3654, 
*0.7188,0.7636,-1.6591,0.7491,0.2532,-0.1620/ 

DATA AD3/-0. 3506, 1.2770, -0.3550, -1.5981, -0.0535, -0.1957, 
*0.4947,-0.1734,0.5484,0.4748,-0.2303,0. 3460/ 

DATA AD4/2. 0684, -1.2777, -0.6224, -0.6550, -0.3275, 0.4715, 
*0.0013,-0.5282,0.0014,-0.2151,-0.3507,-0.3409/ 

DATA DFA/-2. 47, -1.89, -1.65, -0.13, -1.02, -2. 69, -3. 57, 

*-3 .45, -4 . 00, -2 .78 , -7 .0 ,0.83,1.08,0.89, 

*2.14,2.69,5.3,4.09,4.90,5.58,5.39,5.38,5.59,5.89,6.03,6.72, 
*6.51,6.44,6.86,7.14,7.04,6.56,6.29,6.21,-4.2,5.13, 
*4.61,4.85,5.05,3.50,3.46,2.67,2.72,1.68, 2 . 0 , -1 . 2 1 , -1 . 20 , 
*-0. 6,-0. 4 2, -0.3, -0.4, -.87,-1. 35, -1.3, -1.65, -4. 18, 

*-1.3, -3. 99, -4 .29,-3.53/ 

DATA DFB/-2. 1,-1. 4, -1.6, -1.6, -1. 1, -0.2, 0.01, -.15, -.21, 

*-.46, -.74, -.61, .96, 0.2, -1.7, -3. 8, -5. 5, -6. 1,-6. 6, -7. 2, 

*-5.6, -5. 0,-6. 8, -9. 9, -12. 2, -14. 0,-14. 8, -13. 5, -11. 6, -9. 9, 

*-8 .0,-6. 5, -4. 6, -3. 2, -0.8, -3. 9, 2. 8, 2. 50, 1.8,1. ,-6.4, 1.5, 
*2.50,3.9, 5.4, 6.8,3. 00, -2. 1,-3. 2, -2 .8,-3 .7, -6. 6, -9. 70, -10. 7, 
*-8.5, -6. 1,-5. 3, -3. 8, -2. 8, -2. 3/ 

DATA DFC/-1. 5, -1.2, -2. 0,-2. 0,-0. 6, -1.7, -1.4, -1.3, -1.2, 

*-l. 2,-1. 1,-0. 9, -0.3, .1,0. 3, 0.50, -.30, -1.9, -0.8, 0.4, 

*-0. 2,-1. 5, -2. 4, -3. 3, -4. 3, -5. 2, -6. 1,-7. 0,-8. 0,-8. 9, -9. 8, 
*-8.9, -8. 0,-7. 1,-6. 2, -5. 3, -4. 4, -5. 5, -2. 0,-1. 5, -1.0, -0.46, 
*0.05, 0.6, 1.1, 1.6, 2. 1,2. 6, 3. 1,3. 1,2. 9, 2. 1,0. 2, -.7, -0.8, 

*-2.7, -2 .8, -1.5, -1.3, -1.4/ 

DATA DFD/- . 4 , - . 7 , - . 6 , .24,. 67, 1.9, 2. 0,2. 0,2. 2, 2. 3, 2. 2, 

*1.8, 5. 2, 4. 2, 4. 1,4. 7, 5. 1,5. 2, 5. 6, 6. 4, 6. 4, 5. 9, 5. 0,4. 7, 4. 4, 4.1, 


*3. 7, 3. 2, 3 .0,2. 7, 1.9, 1.0, . 

22, . 

64,1. 

3, 2. 5, 2. 6, 2 

.3, 1.9, 1.3, 

*.68, .56, .62, .15, .07,-. 2, . 

07, . 

05,5. 

0 , - . 2 , - . 6 , - 

.3, .1, .5,1.3, 

*1.7, 1.2, .1,-0. 01,-. 3/ 
DATA AZM/0 . , 6 . , 12 . , 18 . , 24 

,30 

.,36. 

, 42. ,48. , 54 

. , 60. ,66. ,72 . , 

*7 8., 84., 90., 96. ,102. ,108. 

,114 

. , 120 

, 126. , 132. 

,138. ,144. ,150. 

*156. , 162. , 168. , 174 . , 180. , 

186. 

,192. 

, 198. ,204. , 

210. ,216. ,222. , 

*228., 234., 240., 246., 252., 

258. 

,264. 

, 270. , 276. , 

282. ,288. , 

*294. ,300. ,306. ,312. ,318. , 

324. 

, 330. 

, 336. ,342 . , 

348. ,354. ,360./ 


DTR=0. 01745329 

LI=12 

LI1=LI-1 
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IF ( ISECOD . EQ . 2 ) THEN 

WRITE (6,32) , ... 

32 FORMAT (1H1/1X, 'SITE ERROR POLYNOMIAL /) 

WRITE (6,3) (AO(I) ,1-1, 4) 

3 FORMAT (IX, # AO*' , 4F10 . 4/) 

DO 5 1-1,12 

WRITE (6,2) I , ADI (I) ,AD2(I) ,AD3(I) ,AD4(I) 

2 FORMAT ( IX , I4,4F10.4) 

5 CONTINUE 

ISECOD— 1 
WRITE(6, 33) 

33 FORMAT (1H1) 

ENDIF 

NDF-4 

C 

C.. CORRECTED BEARING ANGLE Y-AO+A1SINX+A2COSX+ . . . +A11SIN6X+A12COS6X 
C 

STERM-0 . 

CTERM-0 . 

DO 20 J—l , 6 

TJ-DTR*FLOAT(J) 

ITJ-(J-l) *2+1 
ICJ— ( J-l) *2+2 

IF ( IDF . EQ . 1 ) THEN 

STERM-AD1 (ITJ) *SIN (TJ*BRG) +STERM 
CTERM-AD1 (ICJ) *COS (TJ*BRG) +CTERM 
ENDIF 

IF (IDF.EQ. 2) THEN 

STERM-AD2 (ITJ) *SIN (TJ*BRG) +STERM 
CTERM-AD2 (ICJ)* COS ( T J * BRG ) +CTERM 
ENDIF 

IF ( IDF . EQ . 3 ) THEN 

STERM-AD3 (ITJ) *SIN (TJ*BRG) +STERM 
CTERM-AD3 (ICJ) *COS ( TJ * BRG ) +CTERM 
ENDIF 

IF (IDF.EQ. 4) THEN 

STERM-AD4 ( ITJ) *SIN (TJ *BRG) +STERM 
CTERM-AD4 (ICJ) *COS (TJ*BRG) +CTERM 
ENDIF 

20 CONTINUE 

SITER2=A0 (IDF) +STERM+CTERM 
C 

IF (IDF . EQ. IDFTST) SITER2~0. 

C WRITE (6, 25) IDF , BRG , SITER2 

25 FORMAT (1X,I2,2X,2 ( F6 . 2 , 2X) ) 

C DO 120 1*1, LI 

C WRITE (6,35) I , DF2 ( I) , DF3 ( I) , DF4 ( I ) , A2M ( I) 

35 FORMAT (IX, I4,4(F6.2,2X)) 

Cl 20 CONTINUE 

C WRITE (6,399) 

399 FORMAT (IX, 'GOING BACK TO DFDATA ' ) 

909 RETURN 
END 

FUNCTION STDE VB (IDF, BRG ) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


C C 
C FUNCTION STDEVB RECEIVES A DIRECTION FINDER ID AND C 
C AN BEARING ANGLE HAVING A KNOWN SITE ERROR WHICH C 
C IS USED TO COMPUTE THE POLYNOMIAL FORM OF THE BEARING C 
C STANDARD DEVIATION AS A FCN . OF ANGLE FOR USE IN FFIX C 
C C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

COMMON/ DFSTUF/DFE1 (128) ,NS1(128) ,DFR10(128) ,NS10(128) , 
*DFR15 ( 128 ) , NS15 (128) ,SQ1(128) ,SQ10(128) ,SQ15(128) 

COMMON/ KOUNTR/KNT1 , KNT2 , KNT3 , KNT4 , KNTHH,MTKNT, IDFTST, ISECOD 
DIMENSION SDA ( 60) ,SDB(60) ,SDC(60) ,SDD(60) ,AZM(61) 

DIMENSION BO ( 4 ) , BD1 ( 12 ) ,BD2(12) , BD3 ( 12 ) , BD4 ( 12 ) 

DATA BO/1. 63 , -3 .33, -1.47 ,1,94/ 

DATA BD1/0. 8788 , -4 . 9730,-0.4359,0.2113,-0.2401,0.2218, 
*0.4777,0.0991,0.8518,0.1006,0.4715,0.4008/ 

DATA BD2/-2. 1636, 1. 6295,4.9684, -2.5906,-0.0242,1.3654, 
*0.7188,0.7636,-1.6591,0.7491,0.2532,-0.1620/ 

DATA BD3/-0. 3 506, 1. 2770,-0.3550,-1.5981,-0.0535,-0. 1957, 
*0.4947,-0.1734,0.5484,0.4748,-0.2303,0.3460/ 

DATA BD4/2. 0684, -1.2777, -0.6224, -0.6550, -0.3275, 0.4715, 
*0.0013,-0.5282,0.0014,-0.2151,-0.3507,-0.3409/ 

DATA SDA/-2 .47,-1.89,-1.65,-0.13,-1.02,-2.69,-3.57, 

*-3.4 5, -4. 00, -2. 78, -7.0 ,0.83,1.08,0.89, 

*2.14,2.69,5.3,4.09,4.90,5.58,5.39,5.38,5.59,5.89,6.03,6.72, 
*6.51,6.44,6.86,7.14,7.04,6.56,6.29,6.21,-4.2,5.13, 
*4.61,4.85,5.05,3.50,3.46,2.67,2.72,1.68, 2 . 0 , -1 . 21 , -1 . 20 , 
*-0.6, -0.42, -0.3 ,-0.4, -.87, -1.3 5, -1.3, -1.65, -4. 18, 

*-1.3, -3. 99, -4 .29,-3.53/ 
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32 
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2 


DATA SDB/-2. 1,-1. 4, -1.6 # -1.6 , -1.1, -0.2, 0.01 ,--15, -. 21, 

*-. 46 , -. 74 , -. 61, .96, 0.2, -1.7, -3.8, -5. 5, -6. 1,-6. 6, -7,2, 

* -5. 6,-5. 0,-6. 8, -9. 9, -12 .2,-14 .0,-14. 8, -13. 5, -11. 6, -9. 9, 

*-8.0, -6. 5, -4. 6, -3. 2 ,-0.8, -3. 9, 2. 8, 2. 50, 1.8, 1. ,-6.4 ,1.5, 

*2.50,3.9,5.4,6.8, 3. 00, -2 .1,-3. 2, -2. 8 ,-3. 7, -6. 6, -9. 70, -10. 7, 
*-8.5, -6. 1,-5. 3, -3. 8, -2. 8, -2. 3/ 

DATA SDC/-1. 5, -1.2, -2. 0,-2. 0,-0. 6, -1.7 ,-1.4, -1.3, -1.2, 

*-1.2, -1.1, -0.9, -0.3, .1,0. 3, 0.50, -.30, -1.9, -0.8, 0.4, 

*-0.2, -1.5, -2. 4, -3. 3, -4. 3, -5. 2, -6. 1,-7. 0,-8. 0,-8. 9, -9. 8, 

*-8.9, -8. 0,-7. 1,-6. 2, -5. 3, -4. 4, -5. 5, -2. 0,-1. 5, -1.0, -0.46, 
*0.05, 0.6, 1.1, 1.6, 2. 1,2. 6, 3. 1,3. 1,2. 9, 2. 1,0. 2, -.7, -0.8, 

*-2.7, -2. 8, -1.5, -1.3, -1.4/ 

DATA SDD/- .4,— .7,— .6, .24, . 67 ,1.9, 2. 0,2. 0,2. 2, 2. 3, 2. 2, 

*1.8, 5. 2, 4. 2, 4. 1,4. 7, 5. 1,5. 2, 5. 6, 6. 4, 6. 4, 5. 9, 5. 0,4. 7, 4. 4, 4.1, 
*3. 7, 3. 2, 3. 0,2. 7, 1.9, 1.0,. 22,. 64, 1.3, 2. 5, 2. 6, 2. 3, 1.9, 1.3, 

*.68, .56, .62, .15, .07, -.2, .07, .05, 5.0, -.2, -.6, -.3, .1, .5,1.3, 

*1.7, 1.2, .1,-0. 01, -.3/ 

DATA AZM/0. , 6 . , 12 . , 18 . , 24 . , 30 . , 36 . , 42 . , 48 . , 54 . , 60 . , 66 . , 72 . , 
*78. ,84. ,90. ,96. ,102 . , 108. ,114. , 120 . , 126 . , 132 . , 138. , 144. ,150. , 
*156., 162. ,168., 174. ,180., 186., 192., 198., 204., 210., 216., 222., 
*228. ,234. ,240. ,246. ,252, ,258. ,264. ,270. , 276 . , 282 . , 288 . , 

*294 .,300. ,306. ,312. ,318. ,3 24. ,330. ,336. ,342. ,348. ,354. ,360./ 
DTR-0 .01745329 
LI* 12 
LI 1=LI- 1 


IF ( ISECOD . EQ . 2 ) THEN 
WRITE(6, 32) 

FORMAT ( 1H1/ IX , 9 SITE ERROR ST. DEV. POLYNOMIAL * / ) 
WRITE (6,3) (BO ( I ) ,1*1,4) 

FORMAT (IX, ' B0= ' , 4F10.4/) 

DO 5 1=1,12 

WRITE (6, 2) I , BD1 ( I ) , BD2 ( I ) , BD3 ( I ) , BD4 ( I ) 

FORMAT ( IX , 14 , 4F10 . 4 ) 


5 CONTINUE 

ISECOD*-l 
WRITE(6, 33) 

33 FORMAT (1H1) 

ENDIF 

NDF**4 

C 

C.. CORRECTED BEARING ANGLE Y=AO+A1SINX+A2COSX+ . . . +A11SIN6X+A12COS6X 

C 

STERM=0. 

CTERM-0. 

DO 20 J=1 , 6 

TJ-DTR* FLOAT (J) 

ITJ=(J-1) *2+1 
ICJ=(J-1) *2+2 

IF ( IDF. EQ. 1 ) THEN 

STERM-BD1 (ITJ) *SIN (TJ*BRG) +STERM 
CTERM-BDl(ICJ) *COS (TJ*BRG) +CTERM 
ENDIF 

I F ( I DF . EQ . 2 ) THEN 

STERM-BD2 ( ITJ) *SIN ( T J * BRG ) +STERM 
CTERM-BD2 (ICJ) * COS ( TJ * BRG ) +CTERM 
ENDIF 

IF ( IDF. EQ. 3 ) THEN 

STERM-BD3 (ITJ) *SIN (TJ*BRG) +STERM 
CTERM-BD3 (ICJ) *COS (TJ*BRG) +CTERM 
ENDIF 

IF ( IDF. EQ. 4) THEN 

STERM-BD4 (ITJ) *SIN (TJ*BRG) +STERM 
CTERM-BD4 (ICJ) *COS (TJ*BRG) +CTERM 
ENDIF 

20 CONTINUE 

STDEVB-BO (IDF) +STERM+CTERM 
C WRITE (6, 25) IDF, BRG , STDEVB 

25 FORMAT (1X,I2,2X,2(F6.2,2X)) 

C DO 120 1=1, LI 

C WRITE (6, 35) I,DF2(I) ,DF3(I) ,DF4(I) ,AZM(I) 

35 FORMAT (IX, 14, 4 (F6.2,2X) ) 

C120 CONTINUE 

909 RETURN 
END 

SUBROUTINE FFIX ( BRGID, TEMPBR, TM) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
CC THIS PROGRAM IS A VECTOR APPROACH TO DF FIXING CC 

CC MODIFIED AT MSFC TO COMPUTE OPTIMAL LLP SOLNS FROM DF DATA CC 

CC AFTER SITE ERRORS REMOVED FROM DATA (SJG/01-05-86) CC 

C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

IMPLICIT DOUBLE PRECISION (A-H , O-Z ) 

DIMENSION TEMPBR ( 4 ) ,LOD(4) , LOS (4) ,LOM(4) ,LAD(4) ,LAM(4) ,LAS(4) 


DIMENSION BRGSV (50) 

DIMENSION ISUB(50) ,IUSE(50) 
CHARACTER*! N , S , E, W, NORS , EORW 
INTEGER KNAME(50) f BRGID(4) 

DIMENSION SIGX (4 ) , TM ( 4 ) 

DOUBLE PRECISION NXSX, NXSY , NXSZ 
INTEGER SEQ1 , SEQ2 


DOUBLE PRECISION INVAR 
INTEGER SAVE 
DIMENSION SINT (32) 

EQUIVALENCE (SINT(l) ,STAZ(1) ) , (NCNT, NBR) , (IERR, IRES) 

DIMENSION XHISQ (30) 

COMMON K , NST , WBRF , INVAR ( 50 ) , COST ( 50 ) , 

1 SING (50) , C0SG(50) ,STAX(50) , STAY (50) ,STAZ(50) ,ISTA(50) , 

2 BRNX ( 50) , BRNY (50) ,BRNZ(50) ,IL(50) , SAVE (50) ,SWT(50) , 

3 NXSX(50) ,NXSY(50) ,NXSZ(50) # 

4 CHISQ(30) ,DTR,WT(50) ,TX,TY,TZ , El , E2 , E3 , 

5 Cll , C22 , C33 , Cl 3 , C12 , C23 

COMMON/ SOLN/SNKA , SSGNL , TLAT , TLON , RADIUS , AREA , TSOLN , 

* I YDSN , I YR , MON , NDAY , NRS , NBR , IRES , NEVT , SMA , SMI , ORIEN 
COMMON/ KOUNTR/KNT1 , KNT2 , KNT3 , KNT4 , KNTHH, MTKNT, IDFTST, ISECOD 
DATA N,S,E,W/'N' , 'S' , 'E' , 'W'/ 

DATA KNAME/1, 2,3,4, 46*0/ 

DATA LOD/86,86,87,87/, 

* L0H/40, 04,26, 52/ , 

* LOS/ 09 ,37,38, 54/ , 

* LAD/34,35,35,34/, 

* LAM/38,23,50,43/, 

* LAS/57,57,15,00/ 

C ****** USING 20 PER CENT TABLES ********* 

DATA XHISQ/1 .642,3.219,4.642,5.989,7* 289 , 8 . 558 , 9 . 803 , 11.030 , 
B12. 242, 13.442,14.631,15.812, 16.985,18.151,19.311,20.465,21.615, 
C22.760,23.900,25.038,26.171,27.301,28.429,29.553,30.675,31.795, 
D32.912, 34.027, 35. 139, 36. 250/ 

C 

IF(NBR.EQ.O) GOTO 869 
IF (NBR. EQ. 2) THEN 

IF(BRGID(1) • EQ. 0. OR. BRGID(2) .EQ.OJTHEN 
WRITE (6,1191) NBR , BRGID , TEMPBR , TM ( 1 ) ,TK(2) 

1191 FORMAT (IX, 'NULL BRGID FOR NBR -' , 513 , 4F12 . 6, 2X, 2 (F12 . 3 , 2X) ) 

GOTO 871 
ENDIF 
ENDIF 

DTR - .01745329 
TX-0. 

TY-0. 

TZ-0. 

El-0. 

E2-0. 

E3-0. 

Cll-0. 

C22-0. 

C33-0 . 

C13-0. 

C12-0. 

C23-0. 

WBRF— 0 . 

IFLSH-0 

DO 101 11-1,30 
101 CHISQ (II ) - XHISQ (II) 

C NST IS THE NUMBER OF STATIONS IN THE TABLE 

C READ (5,1) NST 

C 1 FORMAT (13) 

NST-4 

DO 20 I— 1 , NST 

C READ THE INDIVIDUAL STATION PARAMETERS 

C READ(5,2) KNAME(I) , LAD, LAM, LAS , LOD, LOM, LOS , SIGX (I ) 

2 FORMAT (A2, 614, IX, F5. 1) 

SLAT - LAD (I) + (LAM (I ) +LAS (I)/60.)/60. 

SLONG - LOD(I) + (LOM ( I) +LOS ( I) / 60 . ) /60 . 

SLAT - SLAT* DTR 
C 

C LONGITUDE INTERNAL TO THE PROGRAM IS NEGATIVE FOR WEST AND 

C POSITIVE FOR EAST. DATA SUBMITED AND PRINTED USES THE OPOSITE 

C CONVENTIONS 

C 

SLONG - -SLONG* DTR 
COSG ( I ) — DCOS ( S LONG ) 

SING (I) - DSIN(SLONG) 

COST ( I ) - DCOS (SLAT) 

C COMPUTE THE STATION VECTOR (ST AX, STAY, STAZ) 

STAZ(I) - DSIN (SLAT) 

STAX ( I ) - COST ( I ) *COSG ( I ) 

STAY ( I ) - COST(I) *SING(I) 
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C WRITE(6,324) STAX(I) ,STAY(I) ,STAZ(I) 

C324 FORMAT ( IX, 3 (F12 . 4 , IX) ) 

C 

C # ## # CALCULATE INVERSE STATION VARIANCE HERE 

C 

C. . SIGX IS THE STANDARD DEVIATION IN DEGREES OF BEARINGS FRO M T HIS 
C 'DF ID#'*I. IT HAS LITTLE OR NO IMPACT ON THE SOLN, BUT THE 
C CONFIDENCE RADII ARE DIRECTLY PROPORTIONAL TO THE VALUE OF 

C SIGX. IN ADDITION, SOME DISTANT SOLNS MAY NOT BE COMPUTED 

C WHEN SIGX IS LARGE (SAY 1 DEG WITH SOLN 500 KM FROM DF) . 

SIGX (I) *1 • 5 

C WRITE (6,3) KNAME ( I ) , LAD ( I ) , LAM ( I ) , LAS ( I ) 

C * , LOD(I) , LOM(I) , LOS (I) , SIGX (I) 

C 3 FORMAT (IX, 12,616, F5. 2) 

20 CONTINUE 
108 FORMAT (A3) 

INDSK-10 
900 CONTINUE 

DO 301 IX* 1, 50 
IUSE(IX) *0 
301 ISUB(IX) *0 
C 

C.. NBR IS THE NUMBER OF BEARINGS IN THIS FLASH 

C 

IFLSH-IFLSH+1 
DO 40 J*1 , NBR 

C READ(9 , 235) BRGID, TEMPBR(J) 

235 FORMAT (IX, 111, IX, F6. 2) 

C WRITE (6 , 241) BRGID, TEMPBR(J) 

C24 1 FORMAT ( IX , 1A1 , IX , F10 . 4 ) 

DO 38 1*1, NST 

C WRITE (6 ,646) BRGID, KNAME (I) 

646 FORMAT (IX, 2 (12 , IX) ) 

IF(BRGID(J) . EQ. KNAME (I) ) GO TO 39 
38 CONTINUE 
C 

C STATION NOT IN TABLE SPACE 


IF(BRGID(J) . NE . KNAME ( I ) ) THEN 
869 WRITE (6,873) NBR, BRGID, TM( BRGID (J) ) 

873 FORMAT ('STATION NOT IN TABLE SPACE? NBR* ',12,', ID- ',414, 

* 2X, F12 . 3 ) 

871 IRES— 1 

GO TO 975 
ENDIF 

39 ISTA ( J) —I 

INVAR (J) —1/ (DTR*SIGX (I ) ) **2 
BRG-TEMPBR ( J ) *DTR 
BRGSV ( J ) — TEMPBR ( I ) 

CBRG - DCOS(BRG) 

SBRG - DSIN(BRG) 

C 

C BEARING PLANE NORMAL VECTOR 

C 

C COMPUTE THE BEARING VECTOR AS THE CROSS PRODUCT OF THE BEARING 

C PLANE NORMALS AND THE SITE VECTORS 

C 

BRNX(J) - CBRG *SING (I) -SBRG *SINT(I) *COSG (I) 

BRNY(J) —CBRG *COSG(I) -SBRG *SINT(I) *SING(I) 

BRNZ(J) - SBRG *COST(I) 

C 

NXSX(J) - BRNY (J) *STAZ (I) -BRNZ (J) *STAY (I) 

NXSY(J) - BRNZ(J) *STAX(I)-BRNX(J) *STAZ(I) 

NXSZ(J) - BRNX ( J) *STAY (I) -BRNY (J) *STAX (I) 

40 CONTINUE 
C 
C 

c 

C USE EXHAUSTIVE REJECTION FOR 10 OR LESS BEARINGS 

C USE SEQUENTIAL REJECTION FOR MORE THAN 10 BEARINGS 

IF (NBR .LE.10) GO TO 50 
CALL SEQ 
GO TO 60 
50 K - NBR 
CALL EXH 

60 IF (IRES .NE. 1) GO TO 800 
C NREJ IS THE NUMBER OF REJECTED BEARINGS 

C K IS THE NUMBER OF BEARINGS USED IN THE FIX 

NREJ - NBR -K 

C WRITE (6, 699) TLAT , TLON 

699 FORMAT (' LAT— ' , F10 . 6, ' LON*',F10.6) 

TLON— TLON 

NORS=N 

EORW-E 
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c WRITE (6,12) LTD , LTM, LTS , NORS , LND, LNM, LNS , EORW, SMA, SMI ,ORIEN f AREA, 

C 1 RADIUS 

12 FORMAT(/, 15X, 'FIX' , 37X, 'S-MAJ AXIS S-MIN AXIS ORIEN ELLIPSE 
1AREA',/15X, 'BPE' f 314 , Al, 2X, 314 ,A1, 8X, 2 (F6. 1 , 5X) , F5 . 1 , 4X, F8 . 1 , 

2/ 1 5X, 'EQI VALENT CIRCULAR RADIUS-' , F6 . 1) 

DO 309 IX-1 , K 

309 IUSE ( ISTA ( IL( IX) ) ) -IUSE ( ISTA ( IL( IX) ) )+l 

C WRITE ( 6 f 307 ) 

307 FORMAT ( 4 (/) # 15X, 'BEARING UTILIZATION' ,//, 15X, ' PDDG' , 3X, ' #SUBMITTE 

ID' , 4X, ' fUSED' , 4X, ' #REJECTED' ) 

DO 303 IX-1, 50 
IF(ISUB{IX) .EQ.O) GO TO 303 
NUSE-ISUB (IX) -IUSE ( IX) 

C WRITE (6,302) KNAME ( IX) , ISUB ( IX) , IUSE ( IX) , NUSE 

302 FORMAT (/, 15X, 15 , 6X, 15 , 5X, 15 , 7X, 15) 

303 CONTINUE 
NUSE-K+NREJ 

C WRITE (6,319) NUSE, K, NREJ 

319 FORMAT (/, 15X, 'TOTAL' , 6X, 15 , 5X, 15 , 7X, 15) 

GO TO 950 
800 IRES - 0 

C WRITE (6, 14) IRES , IFLSH 

14 FORMAT (IX, 11,14, 8X, 'NO FIX' ) 

950 CONTINUE 

C GO TO 900 

975 CONTINUE 
RETURN 
END 

SUBROUTINE EXH 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C c 

C SETS UP FOR BEST POINT ESTIMATE C 

c C 

cccccccccccccccccccccccccccccccccccccccccccccccccccc 
c 

IMPLICIT DOUBLE PRECISION (A-H, O-Z) 

DOUBLE PRECISION NXSX , NXSY, NXSZ 
DOUBLE PRECISION INVAR 
INTEGER SAVE 
DIMENSION KLIM ( 10) 

COMMON K,NST,WBRF, INVAR (50) , COST (50) , 

1 SING (50) , COSG (50) , STAX (50) , STAY (50) ,STAZ(50) , ISTA (50) , 

2 BRNX ( 50 ) , BRNY (50) ,BRNZ(50) ,IL(50) ,SAVE(50) ,SWT(50) , 

3 NXSX ( 50) , NXSY (50) , NXSZ (50) , 

4 CHISQ(30) , DTR, WT(50) , TX ,TY ,TZ , El , E2 , E3 , 

5 Cll , C22,C33,C13,C12,C23 

COMMON/ SOLN/SNKA , SSGNL , TLAT , TLON , RADIUS , AREA , TSOLN , 

* I YDSN , I YR , MON , NDAY , NRS , NBR , IRES , NEVT , SMA , SMI , ORIEN 
COMMON/ KOUNTR/KNT1 , KNT2 , KNT3 , KNT4 , KNTHH , MTKNT , IDFTST , ISECOD 
EQUIVALENCE ( NCNT , NBR ) , ( I ERR , IRES ) 

DATA KLIM/3, 3, 3, 3, 3, 3, 4, 4, 5, 5/ 

DO 330 IP-1,4 

C WRITE (6, 31) STAX ( IP) , STAY (IP) , STAZ (IP) 

31 FORMAT ( IX, 3 (F12 . 4 , IX) , ' IN EXH') 

330 CONTINUE 
10 CONTINUE 

DO 15 I - 1 , K 
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10 CONTINUE 

DO 15 I - 1 , K 


IL( I) - I 
15 CONTINUE 
REJ - 1000. 

C 

c IS STATION UNIQUE? 

C 

20 IS - ISTA ( IL( 1) ) 

DO 21 I - 2 , K 

IF ( IS .NE. ISTA(IL(I) ) ) GO TO 25 

21 CONTINUE 
GO TO 40 

C 

c SET INITIAL WEIGHTS AND CALL BPB 

C 

25 CONTINUE 

DO 26 I«1,K 
WT(IL<I))-INVAR(IL(I)) 

26 CONTINUE 
CALL BPE ( 2 ) 

C 

C IS BEARING SET BETTER THAN PREVIOUS BEST? 

C 

800 CONTINUE 

IF (WBRF .GE. REJ) GO TO 40 
C 

C SAVE THIS CASE 

C 

DO 30 1*1, K 
SAVE ( I ) * IL(I) 

SWT (IL(I) ) *WT(IL(I) ) 

30 CONTINUE 
REJ * WBRF 

C. . THIS FORCES FFIX TO GET SOLN IN ONE PASS IF THERE 
C. . ARE THREE OR LESS BEARINGS IN THE FIX. 

C 

IF (K. LE • 2 ) GOTO 60 

C 

C FIND NEXT COMBINATION 

C 

40 I-K 

41 CONTINUE 

IF ( IL( I) .LT. NBR -K+I) GO TO 45 
IF (I .LE. 1) GO TO 50 
1*1-1 
GO TO 41 

45 J - I 

46 IL( I ) - IL(J)+1 

IF (I .GE. K) GO TO 20 
J - I 
I - 1+1 
GO TO 46 
C 

C ALL COMBINATIONS COMPLETE 

C 

50 WBRF - REJ 

IF (WBRF .LT. CHISQ(K-2) ) GO TO 60 
C 

c TRY SMALLER K 

C 

IF (K .LE. KLIM (NBR) ) GO TO 90 
K - K-l 
GO TO 10 
C 

C OUTPUT RESULTS 

C 

60 CONTINUE 

DO 61 1*1, K 
IL( I ) - SAVE ( I ) 

WT(IL(I) ) *SWT(IL(I) ) 

61 CONTINUE 
CALL BPE ( 1) 

IRES-1 
CALL CONFID 
GO TO 999 

90 IRES - 2 
999 RETURN 
END 

SUBROUTINE SEQ 

f C 
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IMPLICIT DOUBLE PRECIS ION(A-H, 0-2) 

DOUBLE PRECISION NXSX, NXSY , NXSZ 
DOUBLE PRECISION INVAR 
INTEGER SAVE 

COMMON K, NST, WBRF, INVAR(50) ,COST(50) , 

1 SING(50) ,COSG(50) , STAX (50) , STAY (50) ,STAZ(50) ,ISTA(50) , 

2 BRNX{50) , BRNY (50) , BRNZ (50) , IL(50) ,SAVE(50) ,SWT(50) , 

3 NXSX (50) , NXSY (50) , NXSZ (50) , 

4 CHISQ (30) ,DTR,WT(50) ,TX,TY,TZ , E1,B2 , E3 , 

5 011,022,033, C13 , C12 , C23 

COMMON/ SOLN/SNKA , SSGNL, TLAT, TLON , RADIUS , AREA , TSOLN , 

* I YDSN , I YR , MON , NDAY , NRS , NBR , IRES , NEVT , SMA , SMI , ORIEN 
COMMON/ KOUNTR/ KNT 1 , KNT2 , KNT3 , KNT4 , KNTHH , MTKNT , IDFTST, ISECOD 
EQUIVALENCE (NCNT, NBR) , (IERR, IRES) 

K - NBR 
DO 15 I-1,K 
IL(I) - I 
15 CONTINUE 
DO 20 I-1,K 
WT (I) - INVAR (I) 

20 CONTINUE 

21 CONTINUE 
CALL BPE (2) 

C 

C TEST RESULTS 

C 

M « K-2 

IF (K .GE. 33) GO TO 25 
IF (WBRF .LT. CHISQ (M) ) GO TO 50 
GO TO 26 

25 RM-M 

XCHI = DSQRT ( 2 . 0*WBRF) -DSQRT ( 2 . 0*RM-1 . ) 

IF (XCHI . LT. 0.842) GO TO 50 

26 CONTINUE 
C 

C NOT ACCEPTABLE - FIND BEARING TO REMOVE 
IF (K .LE. (NBR+1) /2 . ) GO TO 900 
EMAX - 0. 

DO 30 I-l/K 
IB - IL(I ) 

C FIND THE BEARING HAVING THE LARGEST WEIGHTED ERROR TO THE 

C CURRENT BPE ( X ) 

X - WT (IB) * (BRNX(IB) *TX+BRNY (IB) *TY+BRNZ (IB) *TZ) **2 
IF (0 . LT . TX*NXSX ( IB) +TY*NXSY (IB) +TZ*NXSZ (IB) ) GO TO 60 
X-2*WT(IB)-X 
60 CONTINUE 

IF (X .LE. EMAX) GO TO 30 
EMAX « X 

C SAVE THE INDEX OF THE BEARING TO BE REJECTED 

IS - I 

C WRITE (6, 515) IS 

515 FORMAT (IX, 'REJECTED BEARING ID-' 14) 

30 'CONTINUE 
K - K-l 

C DELETE THE REJECTED BEARING FROM THE LIST 
DO 35 I-IS,K 
IL(I) - IL(I+1) 

35 CONTINUE 
C 

C IS STATION UNIQUE 

C 

II - ISTA(IL(1) ) 

DO 40 I -2, K 

IF (ISTA ( IL(I) ) .NE.I1) GO TO 21 

40 CONTINUE 
C 

C NO FIX 

C 

900 IRES - 2 
GO TO 999 
50 CONTINUE 
CALL BPE ( 1 ) 

IRES— 1 
CALL CONFID 
999 RETURN 
END 

SUBROUTINE BPE (ITER) 

C 

IMPLICIT DOUBLE PRECISION (A-H, 0-Z) 

DOUBLE PRECISION NXSX, NXSY, NXSZ 
DOUBLE PRECISION INVAR 
INTEGER SAVE 
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c 

c 

c 


c 

509 

c 

c 

512 

C 

C513 

C 

C 

c 


COMMON K , NST , WBRF , INVAR ( 50 ) ,COST(50) , 

1 SING (50) , COSG (50) ,STAX(50) ,STAY(50) ,STAZ(50) ,ISTA(50) r 

2 BRNX (50) , BRNY (50) ,BRNZ(50) ,IL(50) ,SAVE(50) ,SWT(50) , 

3 NXSX(50) ,NXSY(50) ,NXSZ(50) , 

4 CHISQ{ 30) , DTR, WT (50) ,TX,TY,TZ,E1, E2 ,E3 , 

5 C11,C22,C33,C13,C12,C23 

COMMON/ SOLN/SNKA , SSGNL, TLAT , TIX)N , RADIUS , AREA , TSOUI , 

*IYDSN,IYR, MON, NDAY,NRS,NBR, IRES, NEVT, SKA, SKI, ORIEN 

COMMON/KOUNTR/KNT1 , KNT2 , KNT3 , KNT4 , KNTHH , KTKNT , IDFTST , ISECOD 

EQUIVALENCE (NCNT, NBR) , (IERR, IRES) 

DATA RFN/ 1.0/ 

DO 100 KK-1, ITER 
Cll* 0. 

C22- 0. 

C33- 0. 

C12- 0. 

C13- 0. 

C23* 0. 

DO 6 1*1, K 
J * IL( I ) 

BUILD THE C MATRIC AS THE PRODUCT OF THE 3 BY N MATRIX OF THE 

WEIGHTED BEARING PLANE NORMAL VECTORS (BRNX, BRNY, BRNZ) AND 

ITS TRANSPOSE 

C11-C11+WT(J) *BRNX(J) **2 

C22-C22+WT(J) *BRNY (J) **2 

C33-C33+WT(J) *BRNZ (J ) **2 

C12-C12+WT(J)*BRNX(J) *BRNY (J) 

C13*C13+WT(J) *BRNX(J) *BRNZ(J) 

C23-C23+WT(J) *BRNY (J) *BRNZ (J) 

6 CONTINUE 

CALL EIGEN (Cll, C22 , C33 , C12 , C13 , C23 , El, E2,E3 ,TX,TY,TZ) 

WRITE ( 6 , 509) C11,C22,C33 

FORMAT (IX, 'C11,C22,C33* ' , 3 (2X, F16 . 4) ) 

WRITE (6 , 509) C12,C13,C23 
WRITE (6, 512) TX,TY,TZ,E1 
FORMAT (IX, 'X,Y,Z*',4 (2X,F12.7) ) 

WRITE ( 6 , 513 ) El 
FORMAT ( ' #- 9 ,F12.7) 

NEED ANTIPODE? 


ICT - 0 
DO 10 1*1, K 
J - IL( I) 

IF (0 . LE * TX*NXSX(J)+TY*NXSY(J)+TZ*NXSZ(J) ) ICT - ICT+1 

10 CONTINUE 

C IF ICT IS K ALL BEARINGS ARE FORWARD 

C IF ICT IS 0 ALL BEARINGS ARE BACKWARD 

IF (K .LE. 2*ICT) GO TO 11 
TX - -TX 
TY - -TY 
TZ - -TZ 

11 CONTINUE 

DO 20 I- 1,K 
J - IL( I ) 

C TDOTT^I^THE COSINE OF THE DISTANCE OF THE BPE TO THE SITE 

TDOTS * TX*STAX (L) +TY*STAY (L) +TZ*STAZ (L) 

IF (TDOTS .GT. .99999) TDOTS - .99999 
IF (TDOTS .LT. -.99999) TDOTS - -.99999 

INSERT RANGE WEIGHTING FUNCTION HERE 
NOTE RFN-1.0 GIVES BEST (SMALLEST) ERROR ELLIPSE FOR LLP 
THE COMPLICATED RFN RELATION ACCOUNTS FOR INCREASED VARIANCE 
IN THE BEARINGS DUE TO SKY WAVE EFFECTS WHEN SOLN CLOSE TO 
THE STATIONS (HIGH INCIDENCE ANGLE TO IONOSPHERE)- THIS IS 
AN IMPORTANT CONSIDERATION IN HF DIRECTION FINDING, NOT HERE. 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

C222 
C 
C 
C 
C 
C 
C 
C 
C 
C 

138 


” ,2F12.3) 


RGE - 34.44* DACOS ( TDOTS ) 

WRITE (6, 222) RGE, TDOTS 
FORMAT (IX, 'RGE AND TDOTS ARE 
IF (RGE .LT. 10.) GO TO 137 
IF (RGE .LT. 1.) GO TO 137 
RFN - .285714+RGE*. 0714256 
GO TO 138 

137 RFN * 3 . -RGE* ( . 402-RGE*. 0204) 

MODIFY THE WEIGHTS BY THE RANGE TO THE CURRENT BPE 
SEE STANSFIELD(1947) FOR WEIGHT 


WT ( J) - INVAR (J) / (RFN* *2* (1 . -TDOTS* *2 ) ) 
20 CONTINUE 
100 CONTINUE 
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WBRF - El 

IF (K .EQ. ICT) RETURN 

IF (0 . EQ. ICT) RETURN 

WBRF - 2*K 

RETURN 

END 

SUBROUTINE EIGEN (Dll , D22 , D33 , D12 , D13 , D23 ,E1, E2 , E3 , TX, TY, TZ) 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c c 

C THIS SUBROUTINE COMPUTES THE SMALLEST EIGENVALUE OF THE MATRIX C C 
C THE EIGEN SUBROUTINE USES NEWTON ITERATION C 

C C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c 

IMPLICIT DOUBLE PRECISION (A-H , O-Z ) 

Cll - Dll 
C22 - D22 
C33 - D33 
C12 =■ D12 
C13 * D13 
C23 - D23 
B - C11+C22+C33 

C - Cll* (C22+C33 ) +C22*C33- (C12**2+C23**2+C13**2 ) 

A-C12/C11 

D-Cll* (C22-C12*A) * (C33-C13*C13/CU- ( (C23-C13*A) **2) 

1 / (C22-C12*A) ) 

IF(D.GT.O) GO TO 13 

El-0 

RETURN 

13 CONTINUE 
X-0 

DO 10 1-1,10 
FX»X**3-B*X**2+C*X-D 
FP-3*X**2-2*B*X+C 
XN-X-FX/FP 
IF(XN.EQ.O) GO TO 12 

IF(ABS ( (XN-X)/XN) .LT. .0001) GO TO 12 
X-XN 

10 CONTINUE 
C WRITE {6 ,107) 

107 FORMAT ( ' 10 ITERATIONS INSUFICIENT TO CONVERGE') 

12 El-XN 

CALL EVECTfEl, D1 1 , D22 , D3 3 , D12 , D13 , D2 3 , TX, TY, TZ) 

RETURN 

END 

SUBROUTINE EVECT (El , Cll , C22 , C33 , C12 , C13 , C23 ,TX, TY, TZ) 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c C 

C THIS SUBROUTINE COMPUTES THE EIGENVECTOR ASSOCIATED WITH THE C 

C EIGEN VALUE El C 

C C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C 

IMPLICIT DOUBLE PRECISION (A-H, 0-2) 

D1 « (Cl 1-El) * (C22-E1 ) -C12**2 
IF (ABS(Dl) .LE. .0001) GO TO 10 
X - Cll* (E1-C22)+C12*C23 

Y - (El-Cll) *C23+C13*C12 
Z - D1 

GO TO 50 

10 D2 - (Cll -El) * (C33-E1) -Cl 3**2 
IF (ABS(D2) .LE. .0001) GO TO 20 
X - (E1-C33) *C12+C13*C23 

Y - D2 

Z - (El-Cll) *C23+C13*C12 
GO TO 50 

20 X - (C22-E1) * (C33-E1) -C23**2 

Y - (E1-C33) *C12+C23*C13 
Z - (E1-C22) *C13+C23*C12 

50 TEMP - DSQRT(X**2+Y**2+Z**2 ) 

TX - X/TEMP 
TY - Y/TEMP 
TZ - Z/TEMP 
RETURN 
END 

SUBROUTINE CON FID 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


C c 
C COMPUTES THE CONFIDENCE ELLIPSE FOR GIVEN DF ANGLE ERROR C 
C WHERE THE ANGLE STANDARD DEVIATION IS IN DEGREES C 
C RETURNS ERROR RADIUS AND AREA FOR THIS ELLIPSE. C 
C C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC GRIG if 1 * A* PAGE IS 
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IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION NXSX, NXSY , NXSZ 
DOUBLE PRECISION INVAR 
INTEGER SAVE 
C 

COMMON K , NST , WBRF , INVAR ( 50 ) , COST (50) , 

1 SING (50) f COSG(50) , STAX (50) , STAY (50) r STAZ(50) ,ISTA(50) , 

2 BRNX (50) , BRNY (50) ,BRNZ(50) ,IL(50) ,SAVE(50) ,SWT(50) , 

3 NXSX(50) f NXSY(50) ,NXSZ(50) , 

4 CHISQ (30) , DTR, WT (50) , TX,TY,TZ, El , E2 , E3 , 

5 Cll C22 C33 C13 C12 C23 

COMMON/ SOLN/SNKA \ SSGNL, TLAT , TLON , RADIUS , AREA , TSOLN , 

*IYDSN , I YR , MON , NDAY , NRS , NBR f IRES , NEVT , SMA , SMI , ORIEN 
COMMON / KOUNTR/ KNT 1 , KNT2 , KNT3 , KNT4 , KNTHH, MTKNT , IDFTST , ISECOD 
C FACT - EARTH RADIUS * DSQRT(-2 LN (1-P)) 

C FACT IS IN UNITS OF NAUT. MI.? FOR P-.5, FACT-4055 

C FACT- 7391 FOR P-0.9 

DATA FACT /4055./ 

MISS-10000. 

c 

TLAT - DAS IN ( TZ ) /DTR 
IF (ABS(TX).LE. l.E-6) GO TO 10 
TLON - ATAN (TY/TX) /DTR 
IF (TX .LE. 0.) TLON - TLON + 180. 

IF (TLON .GT. 180.) TLON - TLON - 360. 

GO TO 15 
10 TLON - 90. 

IF (TY .LE. 0.) TLON - -90. 

15 CONTINUE 
Cll - 0. 

C22 - 0. 

C33 - 0. 

C12 - 0. 

C13 - 0. 

C23 * 0. 

DO 30 1-1 , K 
J - ISTA(IL(I) ) 

C COMPUTE THE BEARING PLANE NORMAL VECTORS (BX, BY, BZ) TO THE BPE 

BX - STAY ( J) *TZ-STAZ ( J) *TY 
BY - STAZ ( J) *TX -STAX ( J) *TZ 

BZ - STAX ( J) *TY - STAY (J) *TX 

C WRITE (6 , 222) COF, BX, BY, BZ ,J, STAX (J) , STAY(J) , STAZ (J) ,TX,TY,TZ 

C222 FORMAT (/IX, 'COF ETC ' , 4F8 . 4 , 12 , 2X, 6 (F8 . 4 , 2X) ) 

BSQAR-BX* * 2+BY* *2+BZ * * 2 
IF (BSQAR. EQ. 0 . ) THEN 
IRES— 0 
GOTO 250 
ENDIF 

COF - WT ( IL ( I ) ) /BSQAR 
C COMPUTE THE WEIGHTED C MATRIX 

Cll - Cll+COF*BX**2 
C22 - C22+COF*BY**2 
C33 - C33+COF*BZ**2 
C12 - C12+COF*BX*BY 
C13 - C13 + COF*BX*BZ 
C23 - C23 + COF*BY*BZ 
30 CONTINUE 

B - (C11+C22+C33J/2. 

C - Cll* (C22+C33 ) +C22*C33- (C12**2+C13 **2+C23**2 ) 

D - DSQRT ( B**2-C) 

E2 - B-D 
E3 - B+D 

C COMPUTE THE CONFIDENCE REGION PARAMETERS 

SMA - FACT/ DSQRT (E2) 

SMI = FACT/ DSQRT ( E3 ) 

CALL EVECT(E2,C11,C22,C33,C12,C13,C23,X,Y,Z) 

IF ( (X*TY-Y*TX) .GT. 0.) Z - -Z 
TEMP - Z/DSQRT( 1. -TZ**2) 

IF (1. .GE. ABS (TEMP) ) GO TO 20 
ORIEN - 0. 

GO TO 22 
20 CONTINUE 

ORIEN - DACOS( TEMP) /DTR 
22 CONTINUE 

AREA - 3. 14159*SMA*SMI 

C CIRCULAR REGION APPROXIMATION (NOT OUTPUT IN THIS VERSION) 

RATIO- (SMI/SMA) 

RADIUS»SMA*( .29294* (RATIO** 2) 063 163*RATIO+. 76996) 

250 RETURN 
END 
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SUBROUTINE NETFIX ( II , AZ1 , 12 , AZ2 , XTL, XTLO) 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


cc 





CC 

cc 

THIS PROGRAM 

COMPUTES A TWO STATION FIX 

CC 

cc 





CC 

cc 

INPUT: 

XSL1 

-LAT OF STATION 1 

(DEGREES) 

CC 

cc 


XSLOl 

-LON OF STATION 1 

(DEGREES) 

CC 

cc 


AZ1 

-OBSERVED BEARING 

(DEGREES) 

CC 

cc 





CC 

cc 


XSL2 

-LAT OF STATION 2 

(DEGREES) 

CC 

cc 


XSL02 

-LON OF STATION 2 

(DEGREES) 

CC 

cc 


AZ2 

-OBSERVED BEARING 

(DEGREES) 

CC 

cc 





CC 

cc 

OUTPUT: 

XTL 

-LAT OF TARGET (DEGREES) 

CC 

cc 


XTLO 

-LON OF TARGET (DEGREES) 

CC 

cc 





cc 

cc 


MODIFIED FOR USE ON EADS 


cc 

cc 


JAN 1987 


cc 


ccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

IMPLICIT DOUBLE PRECISION (A-H, 0-2) 

DOUBLE PRECISION LAD (4) , LAM(4 ) , IAS (4) , L0D(4) , L0H(4 ) , LOS (4) 


DATA LAD/34. 

,35. 

,35. 

,34 

./ 

l , LAM/38 . 

,23. 

,50. 

,43 

•/ 

\ , LAS/57. 

,57. 

,15. 

,00, 

■/ 

, LOD/86 . 

,86. 

,87. 

,87 

./ 

) , LOH/ 40. 

,04. 

,26. 

,52. 


, LOS/ 09 . , 

,37. 

,38. 

,54, 

V 


XSLl-LAD(Il)+(LAM(Il)+(LAS(Il)/60.) )/60. 
XSLOl«LOD(Il)+(LOM(Il)+(LOS(Il)/60.) )/60. 

XSL2-LAD(I2 ) + (LAM(I2 ) + (LAS (I2)/60.))/60. 

XSL02— LOD(I2)+(LOH(I2)+(LOS (I2J/60.) )/60. 

C WRITE (6,901) XSL1 , XSL01 , AZ 1 , XSL2 , XSL02 , AZ2 

C901 FORMAT (IX, 'NETFIX. . . XS LI , XSLOl , AZ1 , XSL2 , XSL02 , AZ2 ' , 6 (F10 . 4 , IX) ) 
PI-3.14159265359 
RAD-57.2957795131 
XSL1-XSL1/RAD 
XSL01-XSL0 1/RAD 
XSL2-XSL2/RAD 
XSL02-XSL02/RAD 

COMPUTE RANGE (RADIANS) BETWEEN SITES 

AZ1-AZ 1/RAD 
AZ2-AZ2/RAD 

WRITE (6,19) XSL1 , XSL01 , XSL2 , XSL02 , AZ 12 , RN 
19 FORMAT( IX, 'BEFORE AZRN ',6F12.4) 

CALL AZRN(XSL1 , XSLOl , X5L2 , XSL02 , AZ12 , RN) 

CALL AZRN (XSL2 , XSL02 , XSL1 , XSLOl, AZ21,RN) 

COMPUTE INCLUDED ANGLES 


ANG1-DABS (AZ12-AZ1) 

IF(ANGl.GT.PI) ANG1— 2 . *PI-ANG1 
ANG2-DABS (AZ21-AZ2 ) 

IF (ANG2 . GT. PI ) ANG2-2 . *PI-ANG2 

COMPUTE 1/2 SUM AND DIFFERENCE ANGLES 

SUM-0 . 5* (ANG2+ANG1 ) 

DIFF-O. 5* (ANG2-ANG1) 

SET UP FOR NAPIER ANALOGY SOLUTION 

SINS— DS IN (SUM) 

SIND-DSIN(DIFF) 

COSS-DCOS (SUM) 

COSD— DCOS (DIFF) 

TANR-DTAN ( 0 . 5*RN) 

COMPUTE LENGTH OF SIDE OPPOSITE STATION 2 


AIN— DAT AN (TANR*SIND/SINS) + DAT AN ( TANR* COSD/ COS S) 

COMPUTE LOCATION OF TARGET 

CALL L0CL0(XSLl,XSLOl,ALN,AZl, XTL, XTLO) 

XTL— XTL*RAD 
XTLO— XTLO ♦RAD 
C WRITE (6, 46) XTL, XTLO 

46 FORMAT (IX, 'NETFIX. .XTL, XTLO 9 , F12 . 7 , 2X, F12 . 7) 
RETURN 
END 
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SUBROUTINE LOCLQ (XLA , XLOT , RT , AT , XLATD, XLOTD) 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


cc ^C 

CC LOCLO COMPUTES THE LATITUDE AND CC 

CC LONGITUDE OF A POINT GIVEN: CC 

CC XLA -LAT OF REFERENCE PT (RAD CC 

CC XLO -LONG OF REF PT(RAD) CC 

CC RT — GCB FROM REF PT TO TARGET CC 

CC AT -OBSERVED AZIMUTH (RAD) CC 

CC OUTPUT: XLATD -LAT OF PT IN QUESTION (RAD) CC 

CC XLOTD -LON OF PT IN QUESTION (RAD) CC 

cc cc 


ccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

IMPLICIT DOUBLE PRECISION (A-H, O-Z) 

PI-3.1415926535 

XAT-AT 

XLO— XLOT 

IF (AT. GT. PI) XAT— DABS (2. *PI-AT) 

XLARL- (PI/2 . ) -XLA 

XLATL-DACOS ( DCOS (XLARL) +DCOS (RT) +DS IN (XLARL) ♦DSIN (RT) *DCOS (XAT) ) 
IF ( XLATL . LT . 0 . ) XLATL-PI+XLATL 
XLATD- ( PI/2 . ) -XLATL 

BETA1-ATAN2 ( (DSIN(0. 5* (RT- XLARL) ) *DCOS ( . 5*XAT) ) , (DSIN( .5* 

C (RT+XLARL) ) +DSIN ( . 5+XAT) ) ) 

BETA2-ATAN2 ( (DCOS ( . 5* (RT-XLARL) ) *DCOS ( . 5*XAT) ) , (DCOS ( . 5* 

C (RT+XLARL) ) *DSIN ( . 5*XAT) ) ) 

GAM— BETA 1+ BET A2 

IF (GAM. LT. 0.0) GAM- (2 . *PI) +GAM 
XLOTD-X LO+G AM 
IF (AT.GE . PI) XLOTD- XLO -GAM 
IF (XLOTD. LT. -PI) XLOTD- (2 . *PI) -XLOTD 
IF (XLOTD. GT. PI) XLOTD-XLOTD- ( 2 . *PI) 

ABVAL— DABS (XLOTD) +DABS (XLO) 

IF ( ( (XLOTD*XLO) .LT.O. ) .AND. (ABVAL. GT. PI) ) XLOTD— 2 . * PI -XLOTD 

IF (XLOTD. GT. 2 . *PI) XLOTD-2 . *PI-XLOTD 

RETURN 

END 

SUBROUTINE AZRN (XSL, XSLO, XTL, XTLO, AZ , R) 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


CC 




cc 

CC 

THIS ROUTINE COMPUTES THE RANGE AND AZIMUTH 

CC 

CC 

FROM ONE 

GEOGRAPHICAL COORDINATE TO ANOTHER 

cc 

cc 




cc 

cc 

INPUT: 

XSL 

-LAT OF POINT S (RADIANS) 

cc 

cc 


XSLO 

-LON OF POINT S (RADIANS) 

cc 

cc 


XTL 

-LAT OF POINT T (RADIANS) 

cc 

cc 


XTLO 

-LON OF POINT T (RADIANS) 

cc 

cc 




cc 

cc 

OUTPUT: 

R 

-RANGE (RADIANS) 

cc 

cc 


AZ 

-AZIMUTH (RADIANS) 

cc 

cc 




cc 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

IMPLICIT DOUBLE PRECIS ION (A-H, O-Z) 

C 

C INITIALIZE EAST WEST FLAG 
C 

C WRITE (6, 33) XSL, XSLO, XTL, XTLO, AZ,R 

33 FORMAT (IX, 'AZRN INPUT.. ',6F12.6) 

IEAWE-0 

PI-3.1415926535 

C 

C A : POLAR ANGLE OF SITE AND TARGET 
C 

A-DABS (XSLO- XTLO) 

C 

C BRANCH IF A IS GREATER THAN PI RADIANS 
IF (A. LE. PI) GO TO 10 
A-2 . *PI— A 
IEAWE-1 
C 

C DS - DISTANCE (RAD) FROM SITE TO NORTH POLE 
C 

10 DS— PI/2 • -XSL 

C 

C DT-D I STANCE (RAD) FROM TARGET TO NORTH POLE 
C 

DT— PI/2 . -XTL 
C 

C AZ1 - . 5* (AZ-B) 

C 



AZ1-ATAN2(DC0S(.5*A) *DSIN ( . 5* (DT-DS) ) f DSIN ( . 5*A) *DSIN ( . 5* (DT+DS) ) ) 


C 

C AZ2 - . 5* (AZ+B) 

AZ2-ATAN2 (DCOS( .5*A) * DCOS ( . 5 * ( DT-DS ) ) ,DSIN( .5*A) *DCOS ( . 5* (DT+DS) ) ) 
AZ-AZ1+AZ2 
C 

C USE LAW OF COSINES FOR SIDES 
C 

R«DACOS(DCOS(DS) +DCOS (DT) +DSIN (DS) +DSIN(DT) +DCOS (A) ) 

C 

C DISTANCE CANNOT BE NEGATIVE 

C 

I F ( R . LT . 0 . ) R— R 
C 

C DETERMINE WHICH WAY TARGET IS FROM SITE 
C 

IF( (IEAWE.EQ. 1 ) .AND. (XSLO. LT. XTLO) ) AZ-2.+PI-AZ 
IF( (IEAWE.EQ. 0) .AND. (XTLO. LT. XSLO) ) AZ-2.+PI-AZ 
C 

C AZ CANNOT BE GT 2*PI OR LT 0. 

C 

IF ( AZ . GE. 2 . *PI) AZ-AZ-2.+PI 
IF (AZ . LT. 0 . ) AZ— AZ 
C 

C CONVERT RADIANS TO KM 
C 

AZ— 2 . *PI-AZ 

RETURN 

END 

SUBROUTINE YDDMY (SYD, DAY, MON, YR) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


C c 

C CONVERT YYDDD (SYD) DA Y , MONTH , YEAR C 

C CONVERT YYDDD TO DAY, MONTH, YEAR (VALID FROM 1901 TO 1999) C 
C DAY - (I) OUTPUT DAY C 

C MON « (I) OUTPUT MONTH C 

C YR - (I) OUTPUT YEAR C 

C c 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

C 

IMPLICIT INTEGER (A-Z) 

DIMENSION DN ( 1 2 , 2 } 

DATA DN/ 3 1,59, 90, 120, 151, 181, 2 12, 243, 273, 304 ,334 ,365, 

* 31,60,91,121,152,182,213,244,274,305,335,366/ 

C 

C 

YY*INT( SYD/ 1000) 

DDD-MOD (SYD, 1000) 

ILY-1 

IF ( MOD ( Y Y , 4 ) . EQ . 0 ) ILY-2 
DO 20 ID-1,12 

IF (DDD. LE. DN (ID, ILY) ) GOTO 21 

20 CONTINUE 

21 MON-ID 
DAY-DDD 

IF (MON . GT. 1 ) DAY-DDD-DN (MON-1, ILY) 

YR-YY 

C WRITE (6, 100) SDY , MON , DAY , YR 

100 FORMAT { IX , ' YYDDD ' ,416) 

RETURN 

END 
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APPENDIX C: CLUSTER ANALYSIS SOFTWARE 


The cluster analysis software consists of programs to read the lightning data from 
a disk file, convert the lightning locations from spherical earth coordinates to rectan- 
gular coordinates within a user defined region from a user defined reference point, 
compute an extrapolation vector, compute seed points for input to the K-means algo- 
rithm, and write the results to an output file. 
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CC CLUSTERING ALGORITHM FOR LTG DATA AND FCSTS OF Y/N AHEAD CC 
CC GRID. READ FROM STAT2A WITH X,Y DATA, NOT LDATA CC 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

COMMON/SO LN/NEVT , TSOLN , I YR , MON , NDAY , TLAT , TLON , SNKA , SSCNL, MRS , 

* SMA , SMI , ORIEN , RADIUS , AREA 

COMMON/ SETGRD/NROW , NCOL, IMIDPT, IRADII , INCR, X, Y f AZ 
COKMON/CM1/LGRID(50,50, 40) , LNUM { 100 , 100) , SGRID(100 , 100) , 

* SNUMflOO, 100) ,AGRID(100, 100) , ANUM( 100, 100) ,A(B000,2) 
INTEGER SNUM.ANUM 

INTEGER IFRQ , I PRINT, IWT, KK, LDCM, LDSWT, LDX , MAXIT , JCOL,NOBS ,NV, NVAR 
DIMENSION IN2 (5) ,CBRG2 (4) ,XBAR{40) , YBAR(40) 

PARAMETER ( IFRQ-0 , 1 PR I NT- 0 , IWT-0 r KX-6 , MAXIT-10 , JCOL-2 , 

* NOBS -6 2 , NV-2 , NVAR- 2 , LDCM- KK, LDSWT- KK, LDX-NOBS) 

INTEGER IC{NOBS) ,IND(NVAR) ,N0B1,NV1 

REAL CM (KK, NVAR) ,NC(KK) , SWT (KK, NVAR) , WSS (KK) , SMTRIX (NOBS , NV) 
CHARACTER I LABEL ( NOBS ) , ICHR (26) 

EXTERNAL KMEAN , WRIRN r WRRRN 

DATA ICHR/ 'A' , / B / ,'C # , , D # , # E / ,'F', # G','H / , , I / ,'J', # K','L' / 'M # , 
♦'N'/O'/P'/Q'/R'/S'/T'/D'/V'/W'/X'/Y'/Z'/ 

DATA IND/1,2/ 

DATA CM (1,1) , CM(1 , 2) , CM (2,1) , CM (2, 2)/ 10. ,60. , 0. , 0*/ 

DATA CM (3,1) , CM (3,2) , CM (4,1) , CM(4 , 2 ) /20. , 0. , -20. , -10-/ 

DATA CM (5,1) , CM (5,2) , CM (6 , 1) , CM( 6 , 2 ) /-70 . , -20 . , -50. , -30 ./ 

DATA FL2LAT, FL2LON, DTR/0 . 604817 , 1* 515038 ,0.017453/ 

DATA CP2IAT, CP2ION/0 .608223,1.515513/ 

DATA BNAIAT,BNALON/0. 632682, 1.510932/ 

DATA CP4 LAT , CP4 LON/ 0.606008,1.515445/ 

DATA ATIME,BTIME, BTIME2/22 3 000 . , 223500. , 223443 ./ 

DATA ALAT , BLAT , ALON , BLON/34 .0,36.5,86.0,90./ 

C 

C. . ADD SYSTEM MOTION VECTOR HERE 

C. . DT— 0 GIVES OBSERVED DATA IF DX , DY USED THEN 

C.. DT WILL BE 5-KIN TRANSLATION OF DATA AND CLUSTERS , DX, DY ARE 

C. . SYSTEM DISPLACEMENT VECTORS FOR 5 -MIN INTERVAL FOR EACH DT 

C. . DT— 2 IS 10 MIN ETC... 

C 

ISTOP— 1 
TINCR— 5000 . 

KZ-1 
DT— 0 . 

DX— 4 . 4 5 *DT 
DY— 7.95*DT 
DO 5 LC-l.KK 

CM ( LC , 1 ) —CM ( LC , 1)+DX 
CM ( LC , 2 ) -CM ( LC , 2 ) +DY 
5 CONTINUE 
XB— 0 . 

YB— 0 . 

I DAY- 15 
KNUMB-0 
IRADII— 200 
I NCR- 5 

NROW-( IRADII *2) /INCR 

NCOL-KROW 
IMIDPT— NROW/2 

C *** READ IN NEXT RECORD FROM DISK FILE STAT2A . 

C 

10 READ( 10 , 5000 ) NNUKB, TSOLN, TLAT, TLON, X, Y 

C CHECK FOR DAY OF INTEREST 
C 

C IF (NDAY . NE. IDAY) GOTO 10 

C IF (TSOLN. GE.ATIME. AND. TSOLN. LT.BTIME. AND. NDAY. £Q. IDAY) THEN 

IF (TSOLN . GE. BTIME2 ) GOTO 37 
IF (TSOLN. LT.ATIME) GOTO 10 

WRITE ( 6 , 5000 ) NNUMB, TSOLN, TLAT, TLON , X,Y 
5000 FORMAT ( 16 , 2X, 5 (F12 . 5 , IX) ) 

IF (TSOLN. GE.BTIME) GOTO 37 

C IF (TLAT. GE.ALAT.AND.TLAT.LT. BLAT. AND. TLON. GE. ALON. AMD. TLON. 

C * LT . BLON) THEN 
C 

c. . CALL AZRN IF NEED CONVERSION OF LAT/ LON TO (X,Y) FROM REF. POINT 
C 

C CALL AZRN(FL2LAT, FL2LON , TARLT ,TARLO, AZ , R) 

C IF (R . LE . IRADI I ) THEN 

C 

C.. IF FORECASTING SYSTEM MOVEMENT, ADD TO X,Y BY DX, DY 
C 

C X— SIN ( AZ ) *R+DX 

C Y-COS(AZ) *R+DY 

C X— SIN (AZ) *R 
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c Y-COS(AZ)*R 

C IF(Y.LE.IOO) GOTO 10 

X-X+DX 
Y-Y+DY 

SMTRIX ( KNUMB+ 1 , 1 ) -X 

SMTRIX (KNUMB+1, 2) -Y 

XB-XB+X 

YB-YB+Y 

KNUMB-KNUMB+1 

GOTO 10 

60 CONTINUE 

WRITE (6,2211) NOBS , NVAR , LDX , IFRQ , IWT , (IND(J) ,J-1,NVAR) , 

* KK MAX IT 

2211 FORMAT (1H1, IX , 'NOBS ETC. 918/) 

C 

C. . ASSEMBLE 5 MIN GRIDS 
C 

37 DO 40 1-1, KNUMB 

JL-20+INT (SMTRIX (1,1)/ 10. ) 

IL— 20-INT(SMTRIX(I , 2)/10. ) 

LGRID(IL,JL,KZ)«LGRID(IL, JL,KZ)+1 
WRITE (6,2882) I , KNUMB, IL, JL, KZ 

WRITE (6 , 2883 ) (SMTRIX ( I , JJ) ,JJ»1,2) , LGRID (IL, JL, KZ) 

2882 FORMAT (IX, 'I, KNUMB, IL,JL, KZ , SMTRIX, LGRID' , 5 (14 , 2X) ) 

2883 FORMAT ( IX, 2F8 . 2 , 2X, 16) 

40 CONTINUE 

C 

C.. FIND MEAN X,Y OF DATA SET FOR SYSTEM MOTION VECTOR 
C 

XBAR ( KZ ) -XB/KNUMB 
YBAR ( KZ ) -YB/KNUMB 

WRITE ( 6 , 2900 ) ATIME , BTIME, KNUMB, KZ , XBAR (KZ) , YBAR (KZ) ,DX,DY,DT 
DO 80 IJ-1,40 
C JI-40-IJ+1 

WRITE (6,2902) IJ, (LGRID(IJ, JJ , KZ) ,JJ-1,40) 

C WRITE ( 6 , 2902 ) JI , (LGRID( JI , JL) , JL*1 , 40) 

80 CONTINUE 

WRITE (6, 2904) (1,1-1,40) 

2904 FORMAT (3X,40(1X,I2) ) 

C IF(ISTOP.EQ.l) STOP 

2900 FORMAT (1H1/ IX, 'INTERVAL- ' , F10. 2 , 2X, F10 . 2 , ' KNUMB- ',16, 

*' KZ— ' , 12 , ' XBAR SYSTEM— ' , F8 . 2 , ' YBAR SYSTEM- ', F8 . 2/ 

*' DX— ' , F8 . 2 , ' DY— ' , F8 . 2 , ' DT-',F8.2/) 

2902 FORMAT ( IX, 12 , 40 ( IX, 12 ) ) 

WRITE (6,1299) NOBS , XBAR (KZ ) , YBAR (KZ) , DT, DX, DY, ATIME, BTIME 
1299 FORMAT ( IX, 'NOBS FOR SYSTEM-' , 14 , 2X, 'XBAR- ' , F10 . 2 , 2X, ' YBAR- ', 

* F10 . 2/2X, ' DT «',F8.2,2X, 'DX- ' , F8 . 2 , 2X, ' DY- ', F8 . 2/2X, ' INTERVAL' , 
*F12 . 3 , 2X, F12 . 3/1H1/) 

C 

DO 2213 JK— 1 , KK 

WRITE (6,2215) JK, (CM(JK, JJ) , JJ-1, 2) 

2215 FORMAT ( IX, 'SEED ' , I4/1X, 2 (F10 . 2 , 2X) ) 

2213 CONTINUE 

C IF(NOBS.NE.-l) STOP 

C 

CALL KMEAN ( NOBS , JCOL , NVAR , SMTRIX , LDX , I FRQ , IWT , IND , KK , MAXIT , 

* CM, LDCM, SWT, LDSWT, IC,NC,WSS) 

CALL WRRRN ( ' CM ' , KK , NVAR , CM , LDCM , 0) 

CALL WRRRN(' SWT' ,KX, NVAR, SWT, LDSWT, 0) 

CALL WRIRN('IC', l,NOBS,IC, 1, 0) 

CALL WRIRN ( 'NC' , 1,KK,NC, 1, 0) 

CALL WRRRN ( ' WSS ' , 1 , KK , WSS ,1,0) 

DO 2230 I— 1 , NOBS 

I LABEL ( I ) — ICHR ( IC ( I ) ) 

C WRITE (6,2244) I, IC(I) , ICHR(IC (I) ) , ILABEL(I) 

C2244 FORMAT (2X, 14, 2X, ' IC-' , 14 , A1 , 2X, Al) 

2230 CONTINUE 

WRITE (6, 6000) 

6000 FORMAT (1H1/6X, 'I' , 6X, 'X' , 10X, ' Y' , 6X, 'CLUSTER'//) 

DO 4000 I— 1 , NOBS 

WRITE (6, 4005) I, (SMTRIX ( I, J) ,J-1,NV) ,IC(I) , ILABEL(I) 

WRITE ( 11 , 4005 ) I, (SMTRIX (I , J) ,J-1,NV) ,IC(I) , ILABEL(I) 

4005 FORMAT (1X,I6,2X,2F10.2,2X,I4, 2X , Al) 

4000 CONTINUE 

KZ-KZ+1 
XB— 0 . 

YB— 0 . 

ATIME— ATIME+TINCR 
BTIME- BTIME+TINCR 
STOP 

993 STOP 
END 
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CC PROGRAM AZRN 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

CC THIS ROUTINE COMPUTES THE RANGE AND AZIMUTH CC 
CC FROM ONE GEOGRAPHICAL COORDINATE TO ANOTHER CC 
CC cc 


CC 

INPUT: 

XSL 

=LAT OF POINT 

S 

(DEGREES) 

cc 

CC 


XSLO 

-LON OF POINT 

S 

(DEGREES) 

cc 

CC 


XTL 

-LAT OF POINT 

T 

(DEGREES) 

CC 

CC 


XTLO 

-LON OF POINT 

T 

(DEGREES) 

CC 

cc 






CC 

cc 

OUTPUT: 

R 

-RANGE (KM) 



cc 

cc 


AZ 

-AZIMUTH (RADIANS) 

CC 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

SUBROUTINE AZRN ( XSL, XSLO, XTL, XTLO, AZ , R) 
IMPLICIT DOUBLE PRECISION ( A-H , O-Z ) 

CC 

C INITIALIZE EAST WEST FLAG 
C 

IEAWE=0 

PI-3.1415926535 

C 

C A : POLAR ANGLE OF SITE AND TARGET 
C 

A-DABS (XSLO-XTLO) 

C 

C BRANCH IF A IS GREATER THAN PI RADIANS 
IF(A.LE.PI) GO TO 10 
A=2. *PI-A 
IEAWE-1 


C 

C DS * DISTANCE (RAD) FROM SITE TO NORTH POLE 
C 

10 DS-PI/2.-XSL 
C 

C DT-DISTANCE (RAD) FROM TARGET TO NORTH POLE 
C 

DT-PI/2 . -XTL 
C 

C AZ 1 = . 5 * ( AZ-B) 

AZ1=ATAN2 ( DCOS ( . 5+A) *DSIN ( . 5* (DT-DS) ) r DSIN ( . 5*A) + DSIN ( . 5* (DT+DS) ) ) 
C 

AZ2-ATAN2 (DCOS ( . 5*A) *DCOS( . 5* (DT-DS) ) , DSIN( • 5*A) *DCOS ( .5* (DT+DS) ) ) 
AZ=AZ1+AZ2 


C 

C USE LAW OF COSINES FOR SIDES 
C 

R-DACOS (DCOS <DS) * DCOS (DT) +DSIN(DS) *DSIN(DT) *DCOS(A) ) 

C 

C DISTANCE CANNOT BE NEGATIVE 
C 

I F ( R . LT . 0 . ) R--R 

C 

C DETERMINE WHICH WAY TARGET IS FROM SITE 

C 

IF( (IEAWE.EQ. 1 ) . AND. (XSLO. LT . XTLO) ) AZ=»2.*PI-AZ 
IF ( (IEAWE.EQ.O) .AND. (XTLO . LT . XSLO) ) AZ*2.*PI“AZ 
C 

c AZ CANNOT BE GT 2*PI OR LT 0. 


IF { AZ . GE. 2 . *PI ) AZ“AZ-2 . *PI 
IF ( AZ . LT . 0 . ) AZ=«-AZ 
C 

c CONVERT RADIANS TO KM 
C 

AZ=2 . +PI-AZ 
R«R*6370. 

RETURN 

END 


184 


ORIGINAL PAGE IS 
OF POOR QUALITY 



cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

CC SEED ALGORITHM FOR LTG DATA AND FCSTS OF Y/N AHEAD CC 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c 


c 

c 

c 

c 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

COMMON / SOLN/NEVT r TSOLN f I YR , MON , NDAY , TLAT , TLON , SNKA , SSGNL, NRS , 

* SMA , SMI , ORIEN , RADIUS , AREA 

COMMON /SETGRD/NROW , NCOL , IMIDPT, IRADII , INCR,X, Y, AZ 
COMMON/ CM 1/LGRID (50,50, 40 ) , LNUM(100, 100) ' *0° ' 

* * SNUM (100,100) , ACRID ( 100 # 100) , ANUM( 100 , 100) ,A(8000,2) 

INTEGER IF^; IPR3ENT, IWT, KK, LDCM,LDSWT,LDX,MAXIT, JCOL,NOBS, NV,NVAR 
DIMENSION IN2 (5) ,CBRG2(4) ,XBAR(40) , YBAR(40) , LSEED ( 40 , 40 , 40) 
INTEGER IC(8000) ,IND(2) ,N0B1,NV1 CUTOT v/(inftft 

REAL CM (80,2) ,NC(80) , SWT (80,2) ,WSS(80) , SMTRIX(8000, 2) 
EQUIVALENCE { LDCM , KK ) , (LDSWT , KK) , (LDX,NOBS) 

EXTERNAL KMEAN, WRIRN , WRRRN 

DATA IFRq!iPRINT,IWT,MAXIT,JCOL,NV,NVAR/0,0,0,10,2,2,2/ 

DATA CM(1. 1),CM(1.2),CM(2,1),CH(2, 2 Z-152.2,-49. 3,-118.8,-27.2/ 

DATA CM(3 , 1) , CM(3 , 2) ,CM(4,1) , CM(4 , 2) /-89 . 0 , -5 . 3 , -49 .9, 2°-V 
DATA CM (5,1) , CM (5,2) , CM (6,1) ,CM(6,2)/-72.2, 15 * 1 ' ^ 

DATA CM (7 , 1) , CM(7 , 2 ) ,CM(8,1) ,CM(8 , 2)/-23 . 2 , 119 . 5, 29 . 5 , i2 9 . 0/ 

DATA FL2LAT, FL2 LON, DTR/0. 604817 ,1.515038,0.017453/ 

DATA CP2LAT,CP2LON/0. 608223, 1.515513/ 

DATA BN ALAT, BNALON/O. 632682 ,1.510932/ 

DATA CP4LAT,CP4LON/0. 606008, 1.515445/ 

DATA ATIME, BTIME , BTIME2/2 14 500 . ,215000. ,223443./ 

DATA ALAT , BLAT , ALON , BLON/ 34.0,36.5,86.0,90./ 


C 

C. . 

c. . 
c. - 
c. . 
c. . 


ADD SYSTEM MOTION VECTOR HERE 

DT-0 GIVES OBSERVED DATA IF DX,DY USED F 

DT WILL BE 5-HIN TRANSLATION OF DATA AND CLUSTERS, DX^Y ARE 
SYSTEM DISPLACEMENT VECTORS FOR 5 -MIN INTERVAL FOR EACH DT 
DT-2 IS 10 MIN ETC. . . 


C 

ISTOP-1 

TINCR-5000. 

KZ-1 

DT-0. 

DX-4 .45*DT 
DY— 7.95*DT 

WRITE (6,2233) ATIME , BTIME , DT, DX, DY 
2233 FORMAT ( 1H1/ IX, 'TIME INTERVAL PROCESSED 
*3F8 . 2/ ) 

DO 5 LC- 1 , KK 

CM(LC,1)-CM(LC,1)+DX 
CM(LC,2)«CM(LC,2)+DY 
5 CONTINUE 
XB-0. 

YB-0. 

I DAY- 15 
KNUMB-0 
IRADII— 200 
INCR-5 


, 2F10 . 2 , ' DT , DX , DY- ' , 


NROW- ( IRADII*2 ) /INCR 

NCOL-NROW 

IMIDPT* NROW/ 2 

C *** READ IN NEXT RECORD FROM DISK FILE STAT2A__. 

10 READ ( 10 , 5000, END-993) NNUMB, TSOLN, TLAT, TLON, X, Y 

IF (TSOLN. GE.BTIME2) GOTO 37 
IF (TSOLN. LT. ATIME) GOTO 10 

C WRITE (6,5000) NNUMB , TSOLN , TLAT , TLON , X , Y 

5000 FORMAT ( 16 , 2X, 5 ( F12 . 5 , IX) ) 

IF (TSOLN. GE. BTIME) GOTO 37 
X-X+DX 
Y-Y+DY 

SMTRIX ( KNUMB+ 1 , 1 ) -X 

SMTRI X ( KNUMB+ 1,2) — Y 

XB-XB+X 

YB-YB+Y 

KNUMB-KNUMB+1 

GOTO 10 

60 CONTINUE 

37 W WRITE(^2211) NOBS , NVAR, LDX, IFRQ, IWT, (IND(J) ,J— 1,NVAR) , 

* KK MAXIT 

2211 FORMAT (1H1, IX, 'NOBS ETC. 918/) 
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c 

C. . ASSEMBLE 5 MIN GRIDS 
C 

DO 40 1-1, KNUMB 

JL-20+INT(SMTRIX(I,l)/10.) 

IL-20-INT(SMTRIX(I,2)/10.) 

LGRID ( IL, JL, KZ ) -LGRID ( IL, JL, KZ ) +1 
LSEED ( IL, JL, KZ ) -LGRID ( IL, JL, KZ) 

40 CONTINUE 
C 

C. . FIND SET OF SEEDS FROM DENSITY GRID 
C 

DO 90 IL-2,39 
DO 90 JL«2,39 

N1-LSEED(IL-1,JL-1,KZ) 

N2-LSBED(IL-1, JL,KZ) 

N3-LSEED ( IL-1 , JL+1 , KZ ) 

N4 -LSEED (IL, JL- 1 , KZ ) 

N5-LSEED(IL,JL,KZ) 

N6— LSEED (IL, JL+1 , KZ ) 

N7-LSEED ( IL+ 1 , JL- 1 , KZ ) 

N8— LSEED( IL+1 , JL, KZ) 

N9— LSEED (IL+1 , JL+1, KZ) 

M1»LGRID(IL-1, JL-1,KZ) 

M2-LGRID ( IL-1 , JL, KZ) 

M3-LGRID( IL-1 , JL+1 , KZ) 

M4— LGRID(IL, JL-1, KZ) 

M5-N5 

M6-LGRID(IL, JL+1 , KZ) 

M7-LGRID (IL+1, JL-1 , KZ ) 

M8— LGRID( IL+1 , JL, KZ) 

M9-LGRID (IL+1, JL+1, KZ) 

I F ( N 5 . LE . N 1 . OR . N 5 . LE . N 2 . OR . N 5 . LE . N 3 . OR . N 5 . LE . N4 . OR . N 5 . LE . N 6 . 

* OR. N5 . LE. N7 * OR. N5 . LE . N8 . OR.N5 . LE .N9) LSEED ( IL, JL, KZ) -0 

IF (M5.LT. Ml. OR. M5. LT.M2.OR.M5.LT. M3. OR. M5.LT.M4. OR. M5.LT.M6) 

* LSEED(IL, JL, KZ) — 0 
90 CONTINUE 

WRITE (6, 190) ATIME , BTIKE 

190 FORMAT (1H1/ IX, 'TIME INTERVAL' , 2 F10 . 2/2X, ' FIRST GUESS FOR SEEDS'/) 
DO 95 IJ-1,40 

WRITE(6, 2902) IJ, (LSEED(IJ, JJ, KZ) ,JJ-1,40) 

95 CONTINUE 

WRITE (6, 2904) (1,1-1,40) 

C 

C.. CONVERT GRID POINT INDEXES INTO X,Y POINTS FOR SEEDING CLUSTERS 
C 

KK-0 

DO 110 1-1,40 
DO 110 J-1,40 

IF ( LSEED ( I , J , KZ ) . GT . 0 ) THEN 
Kl-KK+1 

CM (K1 , 1) —FLOAT ( J-20) *10.+DX 
CM(K1 , 2) —FLOAT (20-1 ) +10 . +DY 

WRITE (6, 2220) I , J , KZ , Kl, LSEED( I , J, KZ) , (CM(KK+1,JJ) ,JJ-1,2) 
2220 FORMAT ( IX , ' I , J , KZ , Kl ',414,' VALUE-', 14,' SEED(X, Y) — ' , 2F8 . 2) 

KK-KK+1 
ENDIF 

110 CONTINUE 
C 

C. . FIND MEAN X,Y OF DATA SET FOR SYSTEM MOTION VECTOR 
C 

XBAR ( KZ ) -XB/KNUMB 
YEAR (KZ) -YB/KNUMB 

WRITE (6, 2900) ATIME, BTIME, KNUMB, KZ , XBAR (KZ) , YBAR(KZ) , DX,DY,DT 
2900 FORMAT ( 1H1/ IX, 'INTERVAL- ' , F10. 2 , 2X, F10. 2 , ' KNUMB- ',16, 

*' KZ— ' , 12 , ' (XBAR, YBAR) SYSTEM- (', F6 . 1 F6 . 1 ,') '/ 

*' DX— ' , F8 . 2 , ' DY— ' , F8 . 2 , ' DT-',F8.2/) 

C 

DO 80 IJ-1,40 
C JI— 40-IJ+l 

WRITE ( 6 , 2902 ) IJ, (LGRID ( IJ, JJ, KZ) ,JJ-1,40) 

2902 FORMAT (IX, 12, 40 (IX, 12)) 

80 CONTINUE 

WRITE (6 , 2904 ) (1,1-1,40) 

2904 FORMAT ( 3X , 4 0 ( IX ,12) ) 

WRITE (6,1299) NOBS , KK , XBAR (KZ) , YBAR (KZ) , DT , DX , DY , ATIME , BTIME 
1299 FORMAT ( IX, 'NOBS , SEEDS-' , 214 , 2X, 'XBAR- ' , F10 . 2 , 2X, ' YBAR- 

*F10.2,2X, 'DT »',F8.2,2X, 'DX- ' , F8 . 2 , 2X, 'DY- ', F8 . 2/2 X, ' INTERVAL' , 
*F12.3,2X,F12.3/) 

IF ( ISTOP . EQ . 1 ) STOP 993 
993 STOP 
END 
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